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DOCUMENTATION  FOR  A  TWO-LEVEL 
DYNAMIC  THERMODYNAMIC  SEA  ICE  MODEL 


William  D.  Hibler  III 


INTRODUCTION 

In  a  recent  article  Hibler  (1979)  formulated  and  tested 
a  two-level  numerical  model  for  the  simulation  of  sea 
ice  circulation  and  thickness  over  a  seasonal  cycle.  Al¬ 
though  developed  for  long-term  simulations,  this  model 
contains  many  features  which  may  prove  useful  in  sea 
ice  forecasting.  In  addition  the  code  may  well  have 
utility  in  studies  of  ice  covered  lakes. 

The  purpose  of  this  report  is  to  document  the  nu¬ 
merical  procedure  and  computer  code  used  in  this  model. 
To  aid  in  using  this  code,  a  complete  listing  and  a 
printout  of  a  double  precision  test  run  are  attached  as 
Appendices. 

In  this  documentation  frequent  reference  will  be 
made  to  the  article  in  the  Journal  of  Physical  Ocean¬ 
ography  that  describes  the  model  (Hibler  1979  here¬ 
after  referred  to  as  JP079).  In  addition  to  describing 
a  variety  of  simulation  results,  this  article  contains  a 
complete  description  of  the  model  equations  and  a 
discussion  of  the  numerical  scheme  and  finite  difference 
code.  Consequently,  various  equations  in  J  P079  will 
frequently  be  directly  referred  to  by  number  in  the  text 
of  this  report. 


BRIEF  DESCRIPTION  OF  MODEL 
Model  structure 

The  overall  structure  (see  Fig.  1)  essentially  consists 
of  three  main  components.  The  first  is  momentum 
balance  which  includes  air  and  water  stress,  Coriolis 
force,  internal  ice  stress,  inertial  forces  and  ocean 
tilt  (sections  2-b-1  and  2-b-2  of  JP079).  Nonlinear 
boundary  layers  for  both  the  ocean-ice  and  air-ice  sur¬ 
face  traction  are  used.  A  key  component  of  this  mo¬ 
mentum  balance  is  the  force  due  to  internal  ice  stress 
(eqs  5  and  6,  JP079).  This  force  is  based  on  a  consti¬ 


tutive  law  which  relates  the  ice  stress  to  the  strain  rate 
and  ice  strength.  For  this  model,  a  viscous-plastic  con¬ 
stitutive  taw  is  used.  Rigid  plastic  behavior  is  approxi¬ 
mated  in  this  law  by  allowing  the  ice  to  flow  in  a  plas¬ 
tic  manner  for  normal  strain  rates  and  to  creep  in  a  lin¬ 
ear  viscous  manner  for  small  strain  rates  (eq  4-1 1  and 
Fig.  1  and  2  in  JP079).  The  treatment  of  ice  as  a  vis¬ 
cous  plastic  fluid  was  largely  motivated  by  the  desire 
to  avoid  the  complexities  associated  with  elastic  plas¬ 
tic  rheologies  (e.g.  Pritchard  1975)  while  still  retain¬ 
ing  plastic  behavior  under  flow. 

The  second  major  component  of  the  model  consists 
of  continuity  equations  (section  2-b-3,  JP079)  describ¬ 
ing  the  evolution  of  the  ice  thickness  characteristics. 
These  equations  are  essentially  a  simplification  to  two 
levels  of  the  Thorndike  et  al.  (1975)  ice  thickness 
destribution  model.  In  this  particular  code  two  cate¬ 
gories  of  ice  are  assumed  :  thin  ice  (<0.5  m  in  thick¬ 
ness)  and  thick  ice  (>0.5  m).  To  keep  track  of  these 
two  categories,  two  parameters  are  calculated:  the 
mean  ice  thickness  per  unit  area  h  and  the  compact¬ 
ness  A  which  is  defined  as  the  fraction  of  area  covered 
by  thick  ice.  The  continuity  equations  describing  the 
evolution  of  these  parameters  (eq  13-16,  JP079)  also 
include  thermodynamic  terms.  For  testing  the  model 
geographically  invariant  growth  rates  as  a  function  of 
thickness  and  time  of  year  were  used. 

A  final  component  is  an  Ice  strength.  For  this 
two-level  case  the  strength  is  taken  to  depend  linearly 
upon  thickness,  and  exponentially  upon  compactness 
(eq  17,  JP079). 

Numerical  scheme 

The  coupled  nonlinear  equations  of  the  model  are 
treated  as  an  initial  value  problem  using  energy  con¬ 
serving  finite  difference  techniques.  A  fixed  grid  is  used 
so  that  the  set  of  equations  may  be  integrated  as  long 
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Figure  1.  Flow  chart  for  two-level  dynamic  thermodynamic  sea  Ice  model. 


as  desired.  Under  this  initial  value  approach,  values  of 
all  parameters  at  all  grid  points  are  required  to  start 
the  integration,  and  boundary  values  of  velocity  are 
required  thereafter.  The  basic  forcing  fields  for  the 
model  consist  of  geostrophic  wind  fields  and  vertically 
integrated  ocean  currents.  Both  inertial  and  momentum 
advection  terms  are  retained  in  the  dynamical  equations. 


Because  of  the  very  strong  ice  interaction  term,  explicit 
integration  of  the  momentum  equations  would  force 
time  steps  to  be  extremely  small  (of  the  order  of  a  few 
minutes).  The  ice  thickness  equations,  on  the  other 
hand,  can  be  explicitly  integrated  over  time  steps  of 
several  days.  To  avoid  this  severe  time  step  limitation 
due  to  the  ice  interaction,  the  momentum  equations 
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are  integrated  implicitly.  This  implicit  integration 
does,  however,  require  a  relaxation  solution  of  a  set  of 
simultaneous  equations  at  each  time  step.  (Details  on 
the  relaxation  scheme  are  given  in  App.  A,  J  P079.) 

Using  this  implicit  scheme  the  inertial  terms  will 
normally  be  significant  for  short  time  steps  (of  the  order 
of  1  hour  or  less)  but  insignificant  for  long  time  steps. 
Although  not  a  formal  stability  requirement,  it  is  wise 
to  choose  time  steps  that  are  small  compared  to  the 
variability  of  the  ice  forcing  because  the  nonlinear  terms 
are  treated  in  a  semi-implicit  manner  in  the  momentum 
balance.  The  effect  of  different  time  step  magnitudes 
on  the  momentum  balance  is  discussed  in  Appendix 
B  of  JP079.  This  appendix  also  contains  some  examples 
which  may  help  the  user  decide  on  time  step  magnitudes 
for  particular  applications.  Because  of  this  implicit 
treatment  of  the  momentum  equation,  the  only  formal 
stability  requirement  is  a  Courant-Friedrichs-Lewy  con¬ 
dition  for  the  advection  terms  in  the  thickness  equations: 

At  <  Ax  (2(w2+v2)] 

A  key  feature  of  the  numerical  scheme  is  a  staggered 
grid  (Fig.  5,  JP079).  This  allows  both  the  ice  strength 
(i.e.  the  nonlinear  viscosities)  and  the  ice  velocities  to 
vary  in  space.  Also  this  configuration  greatly  aids  in 
constructing  energy-conserving  finite  differences.  The 
spatial  finite  differences  used  in  the  numerical  scheme 
are  described  in  some  detail  in  Appendix  A  of  J  P079. 

To  a  large  degree  this  staggered  grid  is  patterned  after 
that  employed  in  primitive  equation  numerical  ocean 
models  (e.g.  Bryan  1 969). 

In  the  time-marching  scheme  (see  eq  21-24,  JP079) 
ice  thickness  and  velocity  are  temporally  staggered. 

In  particular,  the  thickness  characteristics  are  consid¬ 
ered  to  be  defined  at  a  temporal  location  halfway 
between  the  ice  velocities.  This  scheme  (formally 
called  the  forward-backward  procedure)  efficiently 
integrates  coupled  equations  of  this  kind  (e.g.  Mesinger 
and  Arakawa  1976). 

As  mentioned  above,  initial  conditions  at  all  points 
and  ice  velocities  at  the  boundaries  are  thereafter  re¬ 
quired  to  initiate  the  integration  of  the  system  of 
equations  forward  in  time.  The  most  natural  boundary 
condition  is  to  take  the  ice  velocity  to  be  zero  on  the 
boundaries.  This  can  be  done  either  at  a  land  boundary 
or  at  an  ocean  location  where  there  is  no  ice.  Note  that 
the  boundary  condition  does  not  affect  the  ice  motion 
in  such  circumstances  since  in  the  absence  of  ice  the 
strength  is  zero.  More  generally,  as  long  as  the  ocean 
boundaries  are  removed  from  the  ice  edge,  the  coupled 
nature  of  the  model  will  cause  a  natural  ice  edge  boun¬ 
dary  condition  to  be  created  (i.e.  the  ice  strength  drops 
to  zero  naturally  near  the  edge).  However,  it  is  also 


possible  to  form  an  "open"  boundary  condition  by 
setting  the  strength  equal  to  zero  near  a  boundary.  In 
JP079  this  type  of  open  boundary  condition  is  used  at 
the  Spitsbergen-Greenland  passage  to  allow  natural 
inflow  or  outflow.  In  the  computer  code  three  boun¬ 
dary  meshes  are  used  to  define  arbitrary  boundaries  and 
"open”  boundaries.  Islands  are  also  allowed.  Con¬ 
sequently  by  simply  modifying  these  meshes  highly  ir¬ 
regular  boundaries  may  be  taken  into  account. 

Because  of  the  strong  ice  interaction  (which  in  this 
model  is  dissipative  in  nature)  the  momentum  equations 
are  essentially  parabolic  in  form  and  hence  have  few 
numerical  instability  problems  over  long-term  inte¬ 
grations.  (While  there  are  few  such  numerical  problems 
it  should  be  emphasized  that  the  dissipative  ice  inter¬ 
action  terms  are  highly  nonlinear  and  can  lead  to  un¬ 
stable  flow  fields  in  the  absence  of  water  drag.  Such 
a  feature  is  a  physical  characteristic  of  plastic  flow  and 
not  a  numerical  artifact).  However,  in  principle  it  is 
possible  for  the  ice  interaction  to  be  very  small  even 
though  there  may  be  a  finite  ice  mass.  Under  this  situ¬ 
ation,  the  momentum  advection  terms  could  cause  non¬ 
linear  instabilities.  To  ensure  against  such  situations, 
the  bulk  viscosity  parameter  (and  hence  indirectly  shear 
viscosity  as  well)  is  never  allowed  to  drop  below  4.0x  108 
kg  s'1 ,  a  value  which  negligibly  modifies  the  ice  drift. 

Nonlinear  instabilities  over  long-term  integrations 
can  also  arise  from  the  nonlinear  advection  terms  in 
the  thickness  continuity  equations.  To  avoid  this  prob¬ 
lem,  small  biharmonic  and  harmonic  diffusion  terms 
have  been  added  to  the  continuity  equation.  These 
terms  are  grid-size  dependent  in  the  code  and  are 
automatically  made  smaller  for  finer  mesh  grids.  These 
diffusion  terms  are  treated  in  a  mass-conserving  manner 
at  all  boundaries  (see  App.  A,  JP079).  In  practice  these 
diffusive  terms  make  only  small  contributions  (in  the 
JP079  seasonal  simulation  the  average  diffusive  flux 
was  less  than  3%  of  the  average  advective  flux).  It  also 
appears  that  in  the  fully  coupled  model  the  diffusion 
terms  might  be  dispensed  with.  Some  sensitivity  tests 
of  the  effects  of  diffusion  are  discussed  in  Appendix  C 
of  JP079. 

To  integrate  this  model  over  very  long  time  intervals 
requires  a  thermodynamic  code  specifying  ice  growth 
rates  as  a  function  of  thickness  and  time  of  year.  In  the 
present  code,  spatially  invariant  growth  rates  (applicable 
to  the  Arctic  Basin)  from  Thorndike  et  al.  (1975)  are 
used.  For  interested  readers,  time-independent  computa¬ 
tional  procedures  particularly  useful  for  such  growth 
calculations  have  been  developed  by  Semtner  (1976) 
for  perennial  multi-year  ice,  and  Maykut  (1978)  for 
thin  first-year  ice. 


3 


Table  1.  Schematic  listing  of  model  subroutines. 


Level  I  Level  2 


Level  3 


MAIN  (4880-7810) 


ADVECT  (100-930)—  —  DIFFUS  (1190-1490) 
BNDRY  (940-1180) 


/rtLLu  \  ijw-i  /  jv/ 

nr.  Aw  FELLD1  (1740-1990) 

RELAX  (8470-10720)^  p £LLlp  (2000-2330) 

nXMAXM  (12120-12230) 

FORM  (2500-3540) - PLAST  (8000-8460) 

GROWTH  (41 00-4870)— G ROAT E  (38104090) 
GEO  (35S0-38OO) 

MEAN  (7820-7990) 

FGROWP  (2340-2490) 

XSUM  (12240-12350) _ 


COMPUTER  CODE 

A  complete  test  listing  of  the  computer  code  is  given 
in  Appendix  A.  In  this  test  code  an  idealized  5x6  grid 
is  used.  All  arrays  in  this  program  are  indexed  beginning 
at  the  lower  left-hand  corner  point  (1,1).  Columns  are 
denoted  by  constant  x  values,  and  rows  by  constant 
y  values  [i.e.  the  bottom  row  of  an  array  consists  of 
the  points  (1,1),  (2, 1)  etc.].  In  Appendix  8  double 
precision  test  run  results  of  the  code  are  presented. 

In  discussions  of  the  code  below,  line  numbers  of 
specific  statements  will  be  given  in  parentheses.  It  is 
also  noted  that,  while  in  principle  the  code  is  applicable 
to  rectangular  grid  cells,  it  has  only  been  tested  for 
square  grid  cells  (Dx  =  Dy). 

Overall  structure 

The  code  is  essentially  divided  into  three  levels.  The 
first  level  consists  of  the  main  driving  program  (4880- 
7810)  which  sets  up  the  boundaries,  calls  subroutines 
to  numerically  integrate  the  system  of  equations  forward 
in  time,  and  analyzes  and  stores  results.  (A  description 
of  the  flow  in  the  main  program  is  given  below.)  The 
second  level  consists  of  various  subroutines  to  numerically 
advance  equations  in  time  and  to  define  various  nonlinear 
parameters.  These  subroutines,  in  turn,  call  specialized 
subroutines  at  level  3  to,  for  example,  estimate  growth 
rates  or  determine  nonlinear  viscosities  for  plastic  flow 
on  a  prescribed  yield  surface.  Table  1  gives  the  schematic 
location;  Table  2  briefly  summarizes  the  function  of 
each  subroutine  and  Table  3  gives  a  brief  description  of 
the  major  variables.  In  all  subroutines  and  the  main 
driving  program,  a  parameter  statement  is  used  to  de¬ 
fine  dimension  variables.  Thus,  by  a  simple  editing 
command,  array  sizes  can  be  modified  everywhere  in 
the  code. 

Note  from  the  program  listing  in  Appendix  A  that 


the  thickness  and  compactness  grids  are  one  unit  larger 
in  dimension  in  each  direction  than  the  velocity  grid. 
With  the  staggered  grid  configuration  this  allows  thick¬ 
ness  values  to  be  defined  outside  the  velocity  bound¬ 
aries,  a  procedure  which  is  convenient  for  second  order 
differencing  in  the  thickness  equations.  As  mentioned 
earlier  key  features  of  the  code  are  thickness  and  ve¬ 
locity  masks  which  define  the  boundaries.  These  are 
defined  by  subroutine  BNDRY.  Specifically,  UVM 
is  a  mask  for  the  ice  velocity  (1  at  calculated  velocity 
locations,  0  at  boundaries;  islands  are  allowed),  HEFFM 
is  the  thickness  mask  (1  inside  boundaries,  0  outside), 
and  OUT  is  an  outflow  or  open  boundary  mask  (same 
as  HEFFM  except  that  it  is  zero  at  grid  cells  where  an 
open  boundary  is  required).  For  the  test  case  the  masks 
(defined  in  a  file  called  BNDATA-see  line  1005, 
Appendix  A)  create  the  grid  configuration  shown  in 
Figure  2.  At  outflow  grid  cells  the  ice  strength  is 
set  equal  to  zero  and  the  ice  mass  removed.  (All  the 
mass  removed  is  recorded  in  the  main  program.)  Also, 
to  approximately  preserve  second-order  accuracy  for 
advection  into  or  out  of  an  open  cell,  the  mass  and 
compactness  in  outflow  grid  cells  are  estimated  by 
taking  an  average  from  all  the  grid  cells  adjacent  to 
the  open  boundary  (subroutine  MEAN). 

The  input  fields  to  the  model  consist  of  Geostrophic 
wind,  Geostrophic  Ocean  currents  (from  which  tilt  is 
also  calculated  by  the  code)  and  growth  rates.  Other 
parameters  are  simulated  (except  for  initial  conditions 
and  boundary  values).  Some  description  of  the  prepara¬ 
tion  of  input  fields  from  atmospheric  pressure  data  and 
ocean  currents  is  given  in  JP079,  section  3. 

Appendix  B  gives  the  results  of  a  21 -day  test  run  of 
the  complete  code.  For  documentation  purposes  these 
test  results  are  in  double  precision.  In  practice,  however, 
double  precision  is  not  needed  (exu-U  possibly,  in  the 
relaxation  subroutine,  subroutine  RELAX,  for  very 
small  grid  sizes  (say  <25  km)] . 
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Table  2.  Brief  description  of  model  subroutines. 


Name  and  Location 


Function 


MAIN  (4880-7810) 
AOVECT  (100-930) 
DIFFUS  (1190-1490) 
BNDRY  (940-1180) 
RELAX  (8470-10720) 

FELLD  (1500-1730)  J 
FELLD1  (1740-1990)  > 
FELLIP  (2000-2330)  ) 

XMAXM  (12120-12230) 
FORM  (2500-3540) 

PLAST  (8000-8460) 

GROWTH  (4100-4870) 

GROATE  (3810-4090) 
GEO  (3550.3800) 

MEAN  (7820-7990) 

FGROWP  (2340-2490) 


XSUM  (12240-12350) 


Main  driving  program  for  modtl. 

Performs  advection  by  explicit  time  stepping  (modified  Euler). 

Performs  diffusion  by  explicit  forward  time  steps. 

Sets  up  boundary  masks  (UVM,  HEFFM,  OUT). 

Solves  linearized  momentum  balance  with  spatially  varying  bulk  and  shear  viscosities. 
Uses  over-relaxation  techniques. 

Calculate  various  finite  differences  for  use  in  subroutine  RELAX. 

Finds  location  of  maximum  value  of  a  vector. 

Sets  up  forces,  nonlinear  water  drag  coefficients  and  nonlinear  viscosities  (for  plastic 
flow)  for  use  in  each  time  step. 

Calculates  nonlinear  viscosities  based  on  plastic  flow  (specified  by  an  elliptical  yield 
curve)  for  subroutine  FORM. 

Calculates  changes  in  compactness  A  and  mean  ice  thickness  h  due  to  growth,  and 
redistribution  by  ridging. 

Calculates  growth  rates  for  thin  and  thick  ice  in  each  grid  cell. 

Reads  in  geostrophic  wind  data  at  each  time  step  and  ocean  currents  at  the  beginning 
of  the  model  run. 

Estimates  thickness  and  compactness  at  "open  boundary"  grid  cells  based  on  ad¬ 
jacent  values. 

Defines  growth  rates  vs  thickness.  (In  seasonal  runs,  this  subroutine  is  not  used. 
Instead  a  year-long  matrix  of  growth  rates  is  read  in  at  line  5680  and  daily  values 
taken  from  this  matrix  at  lines  6200-6230.) 

Sums  up  values  of  a  vector. _ 


X 

Figure  2.  Boundary  configuration  for  test  grid.  Circled 
grid  points  have  zero  velocities;  shaded  grid  cells  have 
zero  thickness  and  compactness.  Speckled  grid  cell  is 
an  "outflow" cell. 


Description  of  flow  in  main  driving  program 

To  aid  in  understanding  and  using  this  code  it  is  help¬ 
ful  to  trace  through  the  initialization  and  time-stepping 
procedure  in  the  main  driving  program  (lines  4880- 
7810).  Such  a  procedure  also  gives  some  insight  into 
various  subtleties  of  the  numerics  as  well  as  of  the 
computer  code.  In  this  flow  description,  statements 
will  be  referred  to  by  line  number. 

The  first  part  of  the  main  driving  program  (4880- 
6170)  is  concerned  with  initialization.  Early  in  the 
code  (5100-5220)  constants  are  defined  (1-day  time 
steps,  1 25-km  square  grid,  etc.)  and  various  counters 
initialized  (5230-5300).  At  lines  5310-5520  the 
boundary  masks  are  read  in  (using  subroutine  BNDRY) 
and  various  sums  over  the  boundaries  performed  to 
determine  the  number  of  grid  cells  covered  by  ice, 
etc.  The  program  also  prints  out  all  the  boundary  masks 
to  ensure  that  they  are,  in  fact,  the  desired  boundaries. 
At  lines  5530-5640  the  geostrophic  currents  are  read 
in  along  with  the  first  day’s  geostrophic  wind  data 
(using  subroutine  GEO).  Note  that  a  parameter  WK 
(see  line  5540)  is  used  to  cause  GEO  to  either  read  in 
new  currents  or  to  use  old  ones.  At  line  5650-5690, 
a  set  of  growth  rates  is  read  in  (using  subroutine 
FGROWP).  In  the  main  seasonal  simulation  in  JP079, 
FGROWP  is  ignored  and  a  seasonal  matrix  is  read  off 
the  tape  at  line  5680.  However,  for  ice  forecasting 
where  the  growth  rates  do  not  change  very  much, 
FGROWP  with  a  data  statement  could  have  been  used. 
FGROWP  currently  has  January  growth  rates  from 
Thorndike  et  al.  (1 975)-  Growth  rates  are  entered 


Table  3.  Description  of  major  variables  (symbol  in  parentheses  refers  to  variable  symbol  used  in  JP079;  e.g.  h  refers 
to  ice  thickness,  A  to  compactness,  etc.). 


Symbol 


Description 


Location  in  staggered 
grid  (see  Fig.  5,  JPQ79) 


Arrays 


UICE  (u) 

VICE  (v) 

HEFF  (/>) 

AREA  (A) 

ETA  (t,) 

ZETA  (J-) 

AMASS  (m) 

DRAGS  (Ds) 

DRAGA  (0a+M 
GAIRX  ({/g  (x  com¬ 
ponent)) 

GAIRY  (£/g  (y  com¬ 
ponent)) 

GWATX  [t/w  (x  com¬ 
ponent)) 

GWATY  (t/w  (y  com¬ 
ponent)! 

FORCEX  [rx+(3 P/dx)] 

FORCEY  [ry+(d/7dy)] 

UVM 

HEFFM 

OUT 

HCORR 

UICEC  » 

VICEC  I 

OUT1 

OUT2 

UERR1 

VERRl 

HDIFF 

FHEFF  [f(h/A)\ 

FO  |F(o)l 
STRESS  (a) 

HDF 
COR  (X) 

DAIRN 
(=paca|(y  0 
DWATN 

<=pwCwIUw-ul) 

PRESS  (P) 
zmax  (fmax) 

ZMIN  (rmjn) 

GHEFF 

GAREA 


E12  (exy) 


x  component  of  ice  velocity 

y  component  of  ice  velocity 

mean  ice  thickness  per  unit  area 

compactness  (fraction  of  area  covered  by  thick  ice) 

nonlinear  shear  viscosity 

nonlinear  bulk  viscosity 

ice  mass  per  unit  area 

symmetric  nonlinear  water  drag  component  parameter  plus  inertial 
term 

antisymmetric  nonlinear  water  drag  plus  Coriolis  parameter 
x  component  of  geostrophic  wind 

y  component  of  geostroph  ic  wind 

x  component  of  geostrophic  ocean  current 

y  component  of  geostrophic  ocean  current 

x  component  of  external  force  plus  ice  pressure  gradient 
y  component  of  external  forces  plus  ice  pressure  gradient 
boundary  mask  for  velocity  field 
boundary  mask  for  ice  thickness  field 
boundary  mask  for  thickness  excluding  "open  boundaries" 
amount  of  ice  needed  to  remove  "negative  ice  thickness” 
intermediate  ice  velocities  for  use  in  semi-impljcjt  time  step  of  momen¬ 
tum  balance 

difference  between  HEFFM  and  OUT 
dummy  variable  for  use  in  HMEAN 

used  to  obtain  velocity  differences  between  time  steps 

variable  to  keep  track  of  diffusion 
growth  rate  of  thick  ice 
growth  rate  of  thin  ice 

xx,  yy,  and  xy  components  of  Ice  stress  (mks  units) 
variable  to  aid  in  keeping  track  of  diffusion 
Coriolis  parameter 
nonlinear  air  stress  coefficient 

nonlinear  water  stress  coefficient 

ice  pressure  or  strength 

maximum  allowable  bulk  viscosity  value 

minimum  allowable  bulk  viscosity  value 

average  growth  tendency 

area  change  tendency 

xx  strain  rate  component 

yy  strain  rate  component 

xy  strain  rate  component _ 


corner 

corner 

center 

center 

center 

center 

corner 

corner 

corner 

corner 

corner 

corner 

corner 

corner 

corner 

corner 

center 

center 

center 

corner 

center 

center 

corner 

center 

center 

center 

center 

center 

corner 

corner 

corner 

center 

center 

center 

center 

center 

center 

center 

center 


Vectors 

FGROW  (f)  growth  rate  of  ice  at  O.S-m  interval  [There  are  also  a  variety  of  vectors  center 

in  the  program  which  are  equivalenced  to  arrays.  This  was  done  to 
facilitate  sums  and  maximum  point  location  which  can  be  vectorized 
when  vectors  are  used  on  the  Texas  Instruments  Advanced  Scientif  * 

Computer  (ASC)  for  which  this  code  was  initially  written.) 
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Table  3  (cont  d).  Description  of  major  variables  (symbol  in  parentheses  refers  to  variable  symbol  used  in  JP079,  e.g. 
h  refers  to  ice  thickness,  A  to  compactness,  etc.). 


Symbol 

Description 

Constants 

DELTT  AND  DELTAT 

time  step  length  in  seconds 

DELTAX 

grid  size  in  meters 

DELTAY 

grid  size  in  meters 

ERROR 

accuracy  of  relaxation  solution  in  m/s 

DIFFI  (£>,) 

harmonic  diffusion  constant 

A22 

minimum  allowable  value  of  compactness 

HO  </>0) 

demarcation  thickness  between  thin  and  thick  ice 

RHOICE  (/>,) 

ice  density 

FCOR  (\) 

average  Coriolis  parameter 

RHOAIR  (pa) 

air  density 

SINWIN 

sine  of  air  turning  angle 

COSWIN 

cosine  of  air  turning  angle 

SINWAT 

sine  of  ocean  turning  angle 

COSWAT 

cosine  of  ocean  turning  angle 

GRAV  (g) 

acceleration  of  gravity 

ECCEN  («) 

ratio  of  principal  axes  of  plastic  yield  ellipse 

PSTAR  (P«) 

strength  constant 

THETA 

weighting  for  momentum  time  stepping.  THETA  =  1.0  yields  a  back¬ 
wards  time  step;  THETA  =  0.5  a  Crank-Nicholsen  centered  time  step 
and  0.0  a  forward  explicit  time  step. 

in  units  of  centimeters  per  day  and  converted  to  meters 
of  ice  per  second.  Lines  5700-5930  initialize  the  ice 
thickness  and  velocity.  [Note  that  the  velocity,  thick¬ 
ness  and  compactness  functions  are  defined  at  three 
time  levels,  e.g.  H  (/,  /,  K)  K  -  1 , 2,  3.  Of  these  values 
the  K  =  1  level  is  always  the  present  value  and  the  K 
-  2  level  the  previous  time  step.]  As  part  of  this 
initialization,  line  5810  calls  a  subroutine  to  sum  up 
the  thicknesses  at  level  1  over  the  grid,  excluding  the 
“outflow”  grid  cells.  This  is  accomplished  by  using 
the  vector  HEFF1  which  is  equivalenced  to  the  thickness 
array  HEFF.  (Note  that  the  length  of  HEFF1  is  such 
that  only  the  first  level  of  FIEFF  is  included  in  HEFF1 . 
This  procedure  has  been  followed  on  almost  all  the 
equivalence  vectors.)  In  this  initialization,  the  ice  is 
set  to  100%  compactness  with  a  uniform  thickness  of 
3.2967m  (-  3.0/density  of  ice). 

The  remainder  of  the  initialization  (5940-6160)  is 
concerned  with  estimating  the  initial  ice  velocity  field. 
(For  long-term  integrations  this  step  is  not  critical.) 

To  do  this  a  spatially  constant  velocity  is  chosen,  the 
ice  mass  is  set  equal  to  zero  (effectively  removing  the 
inertial  term  from  the  relaxation  solution),  and  the 
parameter  TFIETA  is  set  equal  to  1 .  This  procedure 
gives  a  backwards  time  step-see  eqs  21  and  22,  JP079- 
so  that  the  relaxation  solution  solves  for  the  entire  ice 
velocity  rather  than  portions  of  it,  as  would  happen  in 
a  centered  time  step,  e.g.  TFIETA  =  0.5. 


Flaving  initialized  the  variables,  the  standard  time¬ 
stepping  procedure  is  begun  at  line  61  70.  After  the 
conclusion  of  a  given  time  step,  the  program  returns 
to  line  6180  (more  precisely  to  FORTRAN  statement 
number  100)  and  cycles  through  again.  To  initiate  the 
integration,  growth  rates  applicable  to  the  day  under 
consideration  are  defined  from  the  growth  rate  matrix 
(line  6220).  Since  the  present  wind  field  is  still  valid, 
a  new  wind  field  is  not  called  in  until  the  end  of  the 
time  step— to  be  mentioned  later.  Since  only  fixed  Jan¬ 
uary  growth  rates  are  used  in  the  test  case,  line  6220 
is  commented.  (Note  that  in  this  particular  FORTRAN 
program  comments  are  denoted  by  asterisks.)  In  lines 
6240-6351  thicknesses  and  compactnesses  at  outflow 
points  are  estimated  using  subroutine  MEAN.  Basically 
this  subroutine  uses  grid  cells  containing  ice  adjacent 
to  the  “open”  boundaries  to  estimate  ice  thicknesses 
for  use  in  advection.  All  the  ice  flowing  into  the  “open” 
boundary  cell  is,  however,  explicitly  accounted  for. 

Part  of  this  accounting  is  done  at  line  6351  where  the 
amount  of  ice  in  the  open  cells  (THEFF1 )  at  the  begin¬ 
ning  of  the  time  step  is  calculated. 

At  lines  6360-6550  the  “predictor”  portion  of  the 
momentum  time  step  is  performed  (see  eq  21 ,  JP079). 

A  feature  of  this  process  is  that  at  lines  6390-6420  ice 
velocities  at  the  third  level  [e.g.  UICE  (I, ),  3)]  and 
"centered”  ice  velocities  [UICEC  (I, ))]  are  set  equal 
to  the  present  ice  velocities  (level  one)  before  calling 
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RELAX.  This  is  done  for  two  reasons:  1)  die  initial 
velocity  "guess”  to  start  the  relaxation  procedure  is 
the  ice  velocity  at  the  third  level  and  2)  the  momentum 
advection  term  is  linearized  by  using  UICEC  (see  eq  21 
and  App.  A,  J  P079).  To  carry  out  the  predictor  pro¬ 
cedure,  the  time  step  is  halved  (line  6510)  and  FORM 
is  called.  FORM  uses  the  present  ice  velocity  to  linearize 
the  momentum  equations  as  discussed  earlier. 

In  lines  6560-6760  the  main  forward  time  step  for 
the  momentum  balance  is  carried  out  (see  eq  22,  )P079). 
The  procedure  here  is  the  same  as  in  the  predictor  step 
except  that  a  full  time  step  is  used  and  the  "predicted” 
velocity  is  employed  to  estimate  the  various  nonlinear 
terms.  This  predictor  corrector  procedure  approximately 
centers  the  nonlinear  terms  while  still  maintaining  a  back¬ 
ward  time  step  for  the  "diffusion”  terms. 

Following  the  momentum  time  step,  the  thickness, 
growth  and  decay  terms  are  explicitly  stepped  forward 
in  time  in  lines  6770-6820.  These  operations  correspond 
to  the  finite  difference  equations  (eq  23-24)  in  JP079. 

The  advection  subroutine  (ADVECT)  handles  the  dy¬ 
namical  portions  of  the  compactness  and  thickness  con¬ 
tinuity  equations,  while  the  growth  subroutine  treats 
the  thermodynamic  terms.  (The  growth  subroutine 
would  be  the  logical  place  to  insert  a  more  complete 
heat  budget  subroutine.)  As  noted  at  line  6820,  it  is 
important  to  call  the  growth  subroutine  after  the  ad¬ 
vection  subroutine.  This  is  because  GROWTH  also  han¬ 
dles  the  mechanical  redistribution  by  insisting  that  A  be 
less  than  or  equal  to  1 .  GROWTH  also  corrects  for  prob¬ 
lems  such  as  negative  ice  thicknesses  which  may  arise 
because  centered  differencing  of  second-order  accuracy 
is  used  in  the  advection  code  (such  negative  values  can 
be  avoided  by  upstream  differencing,  but  such  schemes 
are  excessively  diffusive).  Numerical  studies  of  advection 
(Mahlman  and  Sinclair  1977)  show  that  such  negative 
value  removal  substantially  aids  in  tracer  accuracy  pro¬ 
vided  that  it  is  applied  immediately  as  done  here.  Note 
that  in  these  thickness  time  steps  the  latest  values  of  the 
ice  velocity  are  used.  As  mentioned  earlier,  this  effec¬ 
tively  means  that  the  coupled  thickness-momentum 
equations  are  staggered  in  time  (see  eq  21-24,  JP079). 

Once  the  time  stepping  is  completed,  much  of  the 
remainder  of  the  code  is  devoted  to  checking  various 
sums  over  the  array  to  ensure  conservation  and  to  keep 
track  of  the  various  contributions  to  the  thickness  changes. 
In  lines  6830-6933,  the  outflow  ice  and  various  aver¬ 
age  thicknesses  and  growth  sums  are  calculated.  From 
6940-7230  complete  information  over  the  grid  (at  an 
interval  specified  by  lines  6950-6955)  is  printed  out. 

In  lines  7240-7420  velocity  change  information  is 
printed  out  (over  the  same  interval).  In  lines  7440- 
7510  simulated  data  for  each  time  step  are  written 
to  tape  for  later  analysis  (these  write  statements  are 


given  comments  in  the  test  code).  From  lines  751 1  - 
7516  sum  information  is  printed  out  (every  time 
step).  In  lines  7520-7640  new  winds  are  read  in,  and 
the  wind  data  file  is  rewound  if  a  whole  year  is  over. 
Finally  a  check  is  made  on  the  main  time-step  counter 
to  see  if  the  integration  is  finished.  If  not,  the  program 
cycles  back  to  FORTRAN  line  100.  If  finished,  the 
program  terminates  after  writing  out  some  restart  para¬ 
meters  (line  7790)  for  continuing  the  integrations. 

Some  of  the  other  termination  steps  after  line  7650 
are  at  present  superfluous.  They  were  initially  inserted 
for  an  examination  of  “spin  down”  characteristics  of 
the  model. 

Some  comments  on  the  relaxation  subroutine 

Since  the  relaxation  code  is  the  heart  of  the  dynamic 
equations,  some  brief  comments  on  this  subroutine  are 
added  here.  As  presently  written,  the  relaxation  sub¬ 
routine  is  unnecessarily  complex.  There  are  several 
reasons  for  this,  mostly  based  on  computational  efficiency 
considerations.  For  one,  the  relaxation  code  was  sep¬ 
arately  converted  to  single  precision  for  the  seasonal 
runs  discussed  in  JP079  to  speed  up  the  code.  For  this 
purpose  all  variables  except  those  beginning  in  O  were 
defined  in  single  precision.  As  a  result,  the  input  double 
precision  variables  are  converted  to  single  precision  at 
the  beginning  of  the  code  by  a  variety  of  DO  loops. 

Since  such  loops  are  vectorizable  they  take  negligible 
time  on  the  ASC.  In  addition  the  subroutines  FELLD, 
FELLD1  and  FELLIP  were  converted  to  single  precision. 
While  some  changes  may  well  be  useful,  caution  should 
be  exercised  in  any  tinkering  here.  This  is  because  dou¬ 
ble-precision  and  single  precision  variables  are  some¬ 
times  necessarily  mixed  in  the  code  due  to  certain  ASC 
compiler  peculiarities. 

A  second  reason  for  the  complexity  of  the  relaxation 
code,  also  related  to  computational  efficiency,  is  that 
a  large  portion  of  the  relaxation  sweep  was  vectorized 
by  calculating  separately  (in  the  subroutine  FELLD1 ) 
as  much  of  the  finite  difference  code  as  can  be  performed 
over  the  whole  array  without  replacement.  While  in¬ 
creasing  the  complexity  of  the  code,  this  "separation” 
of  the  finite  differences  more  than  doubles  the  speed 
of  the  relaxation  sweep  on  the  ASC. 

A  third  point  is  that  the  relaxation  code  was  made 
general  with  respect  to  backward  or  forward  time  steps. 

In  particular  by  adjusting  THETA  (input  variable 
OTHETA)  backwards,  centered  (Crank-Nicholson)  or 
forward  time  steps  may  be  used.  To  effect  this  generality, 
stresses  at  the  previous  time  step  are  constructed  by 
calling  subroutine  FELLIP.  In  general  the  backwards 
time  step  (THETA  =  1.0)  is  most  stable  but  as  men¬ 
tioned  in  JP079  it  does  overdamp  iti»  «.*:!  waves.  In 
the  Crank-Nicholson  scheme  (THETA  =  0.5,  see,  for 
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example,  Richtmyer  and  Morton  1967)  this  does  not 
occur.  Consequently  such  a  scheme  is  probably 
preferable  at  smaller  time  steps,  say  <  1  hour.  At  longer 
time  steps,  however,  splitting  problems  have  been  found 
to  occur  with  the  centered  scheme.  These  problems 
arise  because  the  inertial  term  is  the  only  term  pulling 
alternating  velocity  variations  together  and  this  term  is 
usually  small  at  large  time  steps.  Some  brief  discussion 
of  this  is  given  in  section  2-C  of  JP079. 

A  final  point  is  that  after  100  sweeps  of  over-relaxa¬ 
tion,  the  subroutine  returns  to  straight  relaxation.  This 
is  because,  while  normally  faster,  the  over-relaxation 
procedure  can  sometimes  diverge. 

CONCLUDING  REMARKS 

As  presently  configured,  this  sea  ice  model  has  direct 
application  to  long-term  simulations.  However,  it  also 
has  many  features  which  may  well  prove  useful  in  ice 
forecasting.  In  addition,  aspects  of  the  code  may  aid 
in  simulating  plastic  flow  in  geophysical  media  other 
than  sea  ice.  In  both  such  applications,  modification 
of  this  code  will  probably  be  helpful.  In  these  cases, 
features  of  the  code  either  relevant  or  necessary  for 
long  integrations  may  be  less  important.  In  addition, 
certain  characteristics  which  play  no  major  role  in  long¬ 
term  simulations  may  prove  important  in  forecasting. 

A  feature  important  for  the  long-term  studies  but 
less  important  for  forecasting  is  the  presence  of  non¬ 
linear  (and  to  a  lesser  degree  linear)  instabilities.  The 
main  examples  of  these  are  alternating  grid  point 
instabilities  occurring  in  the  thickness  and  compactness 
fields  due  to  the  nonlinear  advection  terms.  To  a  lesser 
extent  similar  problems  may  arise  in  the  momentum 
equation  due  to  the  momentum  advection  terms.  Phys¬ 
ically  these  instabilities  arise  due  to  the  nonlinear  ad¬ 
vection  terms  causing  a  cascade  of  energy  to  higher  wave 
numbers.  Due  to  the  finiteness  of  the  grid,  this  energy 
tends  to  pile  up  at  the  Nyquist  wave  number,  producing 
high  wave  number  instabilities.  Because  of  the  com¬ 
pressible  nature  of  sea  ice  flow,  conventional  energy- 
conserving  techniques  cannot  always  prevent  this.  In 
this  code,  this  energy  cascade  is  largely  removed  by 
biharmonic  diffusion.  This  diffusion  preferentially 
damps  out  the  high  wave  number  energy  while  little 
affecting  the  lower  wave  number.  (In  some  sense  this 
procedure  can  be  thought  of  as  a  "poor  man’s  spectral 
model.”)  In  short-term  integrations,  however,  such 
instabilities  never  have  time  to  build  up.  Consequently, 
in  such  applications  there  appears  to  be  no  problem  in 
removing  all  diffusion  terms  in  the  thickness  equations. 
This  removal  will  also  cause  some  slight  linear  instability 
(since  the  advection  terms  are  not  perfectly  centered  in 


time).  However,  over  short  time  intervals  these  instabil¬ 
ities  are  minor.  For  similar  reasons  these  considerations 
also  apply  to  the  minimum  ice  viscosity  used  in  the 
momentum  equations. 

A  converse  feature,  which  is  more  critical  in  the 
short  term,  is  the  problem  of  initialization.  For  long¬ 
term  simulations,  the  model  can  be  integrated  long 
enough  so  that  it  eventually  “forgets”  the  initial  con¬ 
ditions.  For  short-term  forecasts,  however,  the  initial 
conditions  largely  dominate  the  forecast.  How  to  mini¬ 
mize  the  impact  of  poor  initial  conditions  is  an  area 
needing  considerable  research.  One  approach  may  be 
to  use  observed  forcing  data  to  run  the  model  for  a  few 
days  prior  to  the  forecast.  Thus  the  model  may  generate 
some  of  its  own  initial  conditions.  It  may  also  be  that 
update  data  can  be  continually  assimilated  together 
with  this  "initialization”  run. 

Another  difficulty  results  from  specifying  boundary 
conditions  for  localized  regions.  Persistence  may  be 
useful  here.  In  any  case  a  variety  of  these  problems 
need  numerical  examination. 

Finally,  features  of  perhaps  more  general  interest  are 
some  of  the  plastic  characteristics  of  sea  ice  flow.  In  ice 
models  to  date  many  of  these  characteristics  have  been 
suppressed  due  to  the  damping  nature  of  the  oceanic 
boundary  layer.  In  the  absence  of  this  damping,  however, 
the  plastic  flow  is  inherently  unstable  and  becomes  ill- 
defined  without  inclusion  of  the  inertial  terms.  In  par¬ 
ticular,  in  the  absence  ot  water  drag,  the  system  will 
begin  to  accelerate  once  the  plastic  yield  is  exceeded. 

This  behavior  has  been  verified  directly  with  this  code. 
(Initial  results  suggest,  however,  that  this  particular  phe¬ 
nomenon  is  not  critically  dependent  on  time-step  mag¬ 
nitude.)  Consideration  of  spatially  varying  ocean  coup¬ 
ling  effects  with  the  attendant  buildup  of  pressure  terms 
may  well  further  complicate  these  plastic  flow  fields. 

Also  on  shorter  time  scales,  inertial  oscillation  in  both 
the  ice  and  ocean  may  be  important.  Investigation  of 
such  effects  is  needed.  In  these  investigations,  aspects 
of  the  code  described  here  may  prove  useful. 
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APPENDIX  A:  LISTING  OF  COMPUTER  CODE  OF  TWO-LEVEL  MODEL. 


too  SUBROUTINE  ADVECT 1U1CE1 . VICEl. HEFF.  DIFF1,  LAD) 

110*  SPACER 

120  IMPLICIT  DOUBLE  PRECISION  !A-H. 0-Z> 

130  PARAMETER  NX=5,  NY=6, N3» 1 NX  +  1 >  + ! NY+ 1 > , NX  1 =NX+ 1 

140  So  NY1=NY+1,  NXM1=NX-1,  NYM1=NY— 1 

150  DIMENSION  HEFF (NX  1 ,  NY1 ,  3 ) ,  UICE1 ( NX,  NY,  3 ) , VICEl (NX,  NY,  3 > 

160  S<,  UICE ! NX ,  N Y  )  <  VICE  ( NX ,  NY  > 

170  COMMON/STEP/DELTAT, DEL TAX.  DELTAY 

180  COMMON/ ARRAY /H£f-FM ( NX  1 1  NY1  >.  UVM(NX,  NY) 

185  COMMON / D I FF /HD I FF ( NX  1 . N Y 1 1 

190  *  NOW  DECIDE  IF  BACKWARD  EULER  OR  LEAPFROG 

200  LL=LAD 

210  I F  <  LL  EQ  1 )  GO  TO  100 

220  *  BACKWARD  EULER 

230  DELTT=I)El.TAT 

240  K3=2 

250  K2~2 

260  GO  TO  lOl 

270  *  LEAPFROG 

280  100  DELTT =DELTAT*2.  0 

290  K3— 3 

300  K2=2 

310  101  CONTINUE 

320  *  NOW  REARRANGE  H'S 

330  DO  401  J=1  NY 

340  DO  401  1  =  1.  NX 

350  UICE! I. J)=U1CE1 < I. J, 1 ) 

360  VICEl I, J)=VICE1 ( I, J, 1 > 

370  401  CONTINUE 

380  DO  200  J= 1 , NY 1 

390  DO  200  1=1. NX1 

400  HEFF(  I.  J.  3)=HEFF(  I.  J,  2) 

410  HEFF< I.  J,  2> -HEFF ( I, J.  1) 

415  HD IFF < I < J)=0  0 

420  200  CONTINUE 

430  202  CONTINUE 

440  »  NOW  GO  THROUGH  STANDARD  CONSERVATIVE  ADVEC T 1  ON 
450  DELTX=0ELTT/<4  0*DELTAX> 

460  DELTY»DELTT/<4. 0*DELTAY) 

470  DO  210  J=1 , NYM1 

480  DO  210  1=1. NXM1 

490  HEFF ( I  +  l,  J+l.  1 >=HEFF( I  +  l,  J+I , K3 > -DELTX* ( ( HEFF ( l  +  l,  J+l, 2J+HEFF 

500  M 1*2,  J+I.  2) ) * < UICE< 1*1.  J+l >+UICE(I+l, J) ) -<HEFF< l  +  l . J  +  l, 2 ) +HEFF 

510  MI.  J+l.  2)  )*<UICE(  I,  J+l  >+UICE(  I.  J)  >  >  -DELTY*  (  ( HEFF  (  1+1.  J+l,  2)+HEFF 

520  M I  + 1 ,  J+2,  21 )*< VICE! I,  J+l ) +VICE! I  +  l .  J+l ) )-<HEFF( 1  +  1,  J+l. 2 >  +HEFF 

530  S<(  I  +  l.  J,  2)  )*(VICE<  I.  J)+VICE(  I  +  l.  J)  >  ) 

540  210  CONTINUE 

550  *  NOW  DECIDE  IF  DONE 

560  IF  ILL  EQ.  2)  GO  TO  99 

570  IF  (LL  EQ  3)  GO  TO  89 

580  GO  TO  102 

590  89  CONTINUE 

600  *  NOW  FIX  UP  HII.J.2) 

610  00  88  J-1.NY1 

620  DO  88  I-1.NX1 

630  HEFF! I. J.  21-HEFF ( I,  J,  3) 

640  88  CONTINUE 

650  GO  TO  102 

660  99  CONTINUE 

670  «  NOW  DO  BACKWARD  EULER  CORRECTION 
680  DO  220  J-l.NYl 

690  DO  220  1  =  1,  NX  1 

700  HEFF! I.  J.  3)=HEFF(I,  J.  2) 

710  HEFF! I. J. 2 ) =0.  5* ! HEFF ( I , J,  1 >  +HEFF ( I , J, 2> ) 

720  220  CONTINUE 

730  LL*3 
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750  00  TO  202 

760  102  CONTINUE 

770  *  NOW  DO  DIFFUSION  ON  Hd.J.  K3) 

780  DO  240  KD-1.2 

790  00  TO  (241, 242), HD 

800  241  CONTINUE 

810  CALL  DIFFUS (UICE,  VICE,  HEFF,  DIFFI , DELTT ) 

320  00  TO  243 

830  242  CONTINUE 

840  DIFF2=- <  DELTAX++2 )  /DELTT 

850  CALL  DIFFUS<UICE.  UICE,  HEFF, DIFF2, DELTT) 

860  243  CONTINUE 

870  DO  330  J»1,NY1 

880  DO  330  1  =  1, NX  1 

890  HEFFt I, J,  1 ) » <  HcFF ( I , J,  1 ) +HEFF ( I ,  J, 3) >*HEFFM(I. J> 

395  HDIFF  < I,  J ) “HD IFF ( I ,  J)+HEFFd.  J, 3)*HEFFM< I, J) 

900  330  CONTINUE 

910  240  CONTINUE 

920  RETURN 

930  END 

940  SUBROUTINE  BNDRY 

950  *  SUBROUTINE  SETS  UP  BOUNDARY  HASH 

960  IMPLICIT  DOUBLE  PRECISION  (A-H, 0~Z> 

970  PARAMETER  NX=5,  NY*6,  N3= < NX  +  1 >*<NY+1 ) , NX1=NX+1 

980  NY1«NY+1,  NXhl*NX-l.  NYM1-NY-1,  N4=NX*NY 

990  COMMON/ARRAY/H£FFM(NXl,  NY1 ),  UVM(NX, NY) 

1000  COMMON/OUTFLO/OUT (NX1,  NY1 ) 

1005  0PEN(UNIT=1, FILE="BNDATA" ) 

1010  *  READ  IN  VELOCITY  MASH 

1020  READ (UNI T=l. FMT=9 X <UVM< I, U),  1  =  1 , NX  > , J=1 , NY  > 


9  FORMAT (3001  0) 

READ ( UNIT  =  1 ,  FMT=8) ( <HEFFM( I,  J),  I  =  1 , NX  1 ) ,  J= 1 , NY 1 > 

READ(UNIT  =  1,  FMT=8> ( (OUT< I .  U> ,  1  =  1, NX  1 ) . J=1 , NY1 ) 

8  FORMAT (4201.  0) 

RETURN 

END 

SUBROUTINE  DIFFUS(UICE. VICE. HEFF, DIFFI, DELTT > 

*  SPACER 

IMPLICIT  DOUBLE  PRECISION  <A-H, 0-Z) 

PARAMETER  NX  =  5, NY =6.  N3» < NX+1 ) *( NY+ 1 > ,  NX1-NX  +  1 
&, NY1=NY+1, NXM1=NX-1.  NYM1=NY-1 
DIMENSION  HEFF ( NX  1 . NY1 ,  3),  UICE(NX,  NY),  VICECNX,  NY) 

&,  HEFF  1  ( NX  1 ,  NY  1  ) 

COMMON/STEP /DELTAT,  DELTAX,  DELTAY 
COMMON/ ARRAY/HEFFM(NXl,  NY 1 ) ,  UVM(NX,  NY) 

*  SUBROUTINE  DIFFUSES  HEFF, MULTIPLIES  BY  DELT,  AND  PUTS  RESULTS  IN  HEFF 

*  NOW  ZERO  OUT  HEFF1 

DO  210  J=l, NY 1 
DO  210  1  =  1, NX  1 
HEFF 1(1. J  >  =0.  0 
210  CONTINUE 

*  NOW  DO  DIFFUSION 

DELTXX=DELTT*DIFF1 / ( DELTA X*#2 ) 

DELTYY=DELTT*DIFF1/(DEI.TAY**2> 

DO  220  J=2, NY 
DO  220  1=2, NX 

HEFF 1(1, J)=DELTXX*< ( HEFF  < 1*1,  J,  3>-H£FF( X ,  J,  3) )*HEFFM( 1  +  1,  J) 

&  - (HEFF ( I ,  J, 3) -HEFF ( I - 1 ,  J,  3) >+HEFFM( I  — l,  J) ) 

\  +DELTYY* ( (HEFF  < I ,  J+l ,  3) -HEFF ( I,  J,  3) )*HEFFM( I,  J+l > 

V.  -(HEFF<  1,  J,  3 )  -HEFF  (  I ,  J-l,  3)  )  #HEFFM  (  I,  J-l  )  ) 

220  CONTINUE 

DO  260  J=l, NY1 
DO  260  1=1, NX1 
HEFF ( I . J, 3)=HEFF1 ( I,  J) 


1470 

1480 

1490 

2340 

2350 

2360 

2370 

2380 

2385 

2390 

2400 

2410 

2420 

2430 

2440 

2450 

2460 

2470 

2480 

2490 

2500 

2510 

2520 

2530 

2540 

2550 

2560 

2570 

2580 

2590 

2600 

2610 

2620 

2630 

2640 

2650 

2660 

2670 

2680 

2690 

2700 

2710 

2720 

2730 

2740 

2750 

2760 

2770 

2780 

2790 

2800 

2810 

2820 

2830 

2840 

2850 

2860 

2870 

2880 

2890 

2895 

2900 

2910 

2920 

2930 

2940 


260  CONTINUE 
RETURN 
END 

SUBROUTINE  FQRCtWP 

*  SUBROUTINE  TO  SET  UP  CANONICAL  GROWTH  VALUES 

IMPLICIT  DOUBLE  PRECISION  (A-H, D-Z' 

PARAMETER  NX-5.  NY-t,,  N3-.NX+1  >*(NY*1  ).  Wtl«NX+l 
l<,  NY  I  —NY ♦  1  .  NXM1-NX-1 .  NYM 1  “NV -  1 
DIMENSION  FGCOW ( 21 ) 

COMMON/OPOWTH /TGPOW < 2 1 

DATA  FGGOW/12  O9D+00.  1  95D+00.  .  46D+00.  370*00.  31D+00.  27D+0C- 

«...  21D+00,  14D+00,  09D+0G,  03D+00.  0  0  -  02D+00.  -  03D+00, -  040*00 
&• -. 05D+00  -  05D*00.  -  060*00, -  06D+00.  -  0oD*00,  -  060+00  -  060+00/ 
DO  lOl  1-1,21 

FCROWf I >»FGGOW( I >/(8  64D+06) 

PRINT  1,1, FCROW ( I ) 

101  CONTINUE 
1  FORMAT < IX.  15, G20  12 > 

RETURN 

END 

SUBROUTINE  F0RMCU1CE,  VICE,  ETA,  ZETA,  DRAGS 

DRAGA,  GAIRX,  GAIRY,  GWAT x ,  GWATY,  FORCEX,  FORCEY,  HEFF,  AMASS,  AREA) 
IMPLICIT  DOUBLE  PRECISION  (A-H, 0-Z) 

PARAMETER  NX-5,  NY-6,  N3=(NX*1 )* (NY+1 >,  NXl-NX+1 
Si.  NY  1  —NY*  1 ,  NXM1  “NX  — 1 ,  NYM1  “NY-1,  N4«NX*NY 

*  PROGRAM  FORMS  BASIC  INPUT  PARAMETERS  FOR  RELAXATION 

DIMENSION  UICE ( NX,  NY,  3),  VICECNX,  NY,  3), ETACNX1,  NY I ) 

Si,  ZETACNX1.  NY1  >,  DRAGS!  NX  .  NY),  DRAGA  (NX.  NY),  CAIRX(NX,  NY) 
fc,  GAIRY  (NX,  NY),  GWATX(NX,  NY  ),  GWATY  (  NX,  NY  > ,  HF.FF  ( NX  1 ,  NY1 ,  2; 

Si,  COR  (NX,  NY),  FORCEX  (NX,  NY  ),  FORCEY  ( NX.  NY  > ,  AMASS  ( NX .  NY) 

Si,  DAIRN(NX,  NY) ,  DWATN(NX,  NY ) ,  AREA  <  NX  1 ,  NY  I,  3),  PRESS  (NX  I,  NY  1  ) 

Si,  ZMAX  ( NX  1 ,  NY  1  ) ,  ZMIN(NXI,  NYJ  > 

COMMON/STEP/DELTAT.  DEL TAX,  DELTAY 
COMMON/ARRAY/HEFFM(NXl,  NY1 >,  UVMtNX,  NY) 

COMMON/OUTFLO/OUT  <  NX 1 .  NY 1 ) 

*  FIRST  SET  UP  BASIC  CONSTANTS 

DWAT=0.  59D+-00 
DAIR-O.  01462D+00 
RHOICE-O  91D+03 
FCOR—1 .  46D-04 
RHOAIR—1  3D+00 
SINWIN-0  4226D+00 
COSWIN-O  9063D+00 
SINWAT-O.  4226D+00 
CQSWAT-O.  9063D+00 
GRAV=9.  832D+00 
ECCEN-2  0 

*  20  DEG  GIVES  SIN  EQUAL  TO  342)2 

*  NOW  SET  UP  MASS  PER  UNIT  AREA  AND  CORIOLIS  TERM 

DO  101  J-l.NY 
DO  101  1  =  1,  NX 

AMASS ( I ,  J)-RH0ICE*0  25* ( HEFF ( I , J,  1 ) 

S<  +HEFF  ( I  +■  1 ,  J,  1  )  +HEFF  (  I ,  J*l,  1  )  +HEFF  (  I  +  1 ,  J+i,  l  >  ) 

COR ( I , J ) -AMASS ( I ,  J)*FCOR 
101  CONTINUE 

*  NOW  SET  UP  NON  LINEAR  WIND  AND  WATER  DRAG 

DO  99  J-l.  NY 
DO  99  1=1, NX 

DAIRN( I. U >  — RHOA I H*  0012D*00*DSQRT (GAIRX ( I, J > »*2*GAIRY ( I , J>**2) 
DWATN (  I,  J>«5.  5D*-00*DSQRT(  (UICE(  I.  J.  I  >-GWATX(  I,  J)  )**2 
Si  +•( VICE!  I .  J,  I ) -GWATY ( I .  J )  )**2) 

99  CONTINUE 

*  NOW  SET  UP  SYMMETTR 1C  DRAG 

DO  102  J-l, NY 
DO  102  I-l.NX 

DRAGS ( I,  J) -DWATN ( I .  J > *C OSWAT 
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2950  102  CONTINUE 

2960  *  NOW  SET  UP  ANTI  SYMMETRIC  DRAG  PLUS  CORIOLIS 
2970  DO  103  J*  t ,  NY 

2980  DO  103  I  —  1 .  NX 

2990  DRAGA< I. Jl-DWATN! I.  V  >  *SINWAT +-COR < I . J> 

3000  103  CONTINUE 

3010  *  NOW  SET  UP  PORCINO  FIELD 
3020  DO  104  J-l.NY 

3030  DO  104  I-l.NX 

3040  *  FIRST  DO  WIND 

3050  FORCEX  < I ,  J ) -DAIRN ( I . J>*< COSWIN#CAIRX ( I , J> 

3060  S<-SINWIN*CAIRY! I. J) ) 

3070  FORCEY < I > J) -DAIRN! I,  U> *! S1NWIN*GAIRX 1  I ,  J) 

3080  &  +CQSW1N*0AIRY(  I.  J) ) 

3090  *  NOW  ADD  IN  CURRENT  FORCE 

3100  FORCEX <  I ,  J) -FORCEX!  I.  U>  +DWATN! I,  J )  # ! COSWAT «GWAT X !  I .  J) 

3110  S<  -SINWAT*GWATY<  I,  J)  ) 

3120  FORCEY <  I-  JJ-FORCEY!  t.  J)-rDWATN(  I.  U)*(SINWAT*GWATX<  I.  J) 

3130  S.  ♦COSWAT*OWATY<I,  J)  ) 

3140  *  NOW  ADD  IN  TILT 

3150  FORCEX!  I.  U) -FORCEX  <  I.  J>-COR<  I.  J>*OWATY(  I,  J) 

3160  FORCEY! 1/  J)=FORCEY! I.  J)+COR( I.  J>*GWATX< I,  J) 

3170  104  CONTINUE 

3180  *  NOW  SET  UP  ICE  PRESSURE  AND  VISCOSITIES 
3190  *  FIRST  SET  UP  CONSTANTS 
3200  PSTAR-5  0D+03 

3210  EPSIL-5  787D-03 

3220  #  3  787D-08  IS  0  3  PER  CENT  PER  DAY  STRAIN  RATE 

3230  ZETAC-0  IS 

3240  ETAC-ZETAC/4  0 

3250  *  NOW  SET  UP  VALUES 

3260  DO  105  J— 1 1 NY  1 

3270  DO  105  1*1, NX1 

3280  PRESS!  I,  J > =PSTAR*HEFF !  I,  J.  I  >*DEXP!-20  0*<  I.  O-AREA I  I ,  J,  1  )  )  ) 

3290  ZMAX!I.  J)*!S  OD+12/!2.  OD+04 > > *PRESS ! I . J) 

3300  ZMIN< I, J>-4.  0D+08 

3310  *  ZETA! I,  J)=ZETAC*PRESS< I- Ul/EPSIL 

3320  *  ETA! I .  J ) =ETAC*PR£SS ! I ,  J) /EPSIL 
3330  *  ZETA! I,  J) “OMAX 1 ! 4.  00+08.  ZETA! I >  J ) ) 

3340  *  ETA! I ,  J ) -DMAX 1 ! 1 .  0D+08.  ETA! I.  U) ) 

33S0  105  CONTINUE 

3360  CALL  PLAST!UICE.  VICE,  PRESS.  ETA,  ZETA, ECCEN,  ZMAX,  ZMIN) 

3370  *  NOW  SET  VISCOSITIES  AND  PRESSURE  EQUAL  TO  ZERO  AT  OUTFLOW  PTS 
3380  DO  106  U-1.NY1 

3390  DO  106  1*1, NX  1 

3400  PRESS! I.  J) -PRESS! I,  J)*OUT! I,  J) 

3410  ETA! I,  J)*ETA< I, U)*OUT< I,  U> 

3420  ZETA! I. J  >  *ZETA< I ,  J)*OUT! I. J) 

3430  106  CONTINUE 

3440  *  NOW  CALCULATE  PRESSURE  FORCE  AND  ADD  TO  EXTERNAL  FORCE 
3450  DO  107  J-l.NY 

3460  DO  107  1=1, NX 

3470  FORCEX! I, J) -FORCEX! I. J)-!0.  25/DELTAX ) * 

3480  Si  !  PRESS !  I  + 1 .  V  >  +PRESS !  I  ♦  1 ,  U-M) -PRESS!  I,  J)  -PRESS !  t  >  U+ 1 )  ) 

3490  FORCEY!!. J)»FORCEY! I. J>-!0  25/DELTAY > * 

3500  Si  ! PRESS!  I,  J«-l  IMPRESS!  1  +  1,  U+l  ) -PRESS!  I,  U) -PRESS!  1  +  1,  J) ) 

3510  #  NOW  PUT  IN  MINIMAL  MASS  FOR  TIME  STEPPING  CALCULATIONS 
3520  107  CONTINUE 

3530  RETURN 

3540  END 

3550  SUBROUTINE  GEQ!GAIRX, GAIRY, GWATX, GWATY, WK) 

3560  *  IF  WK  LT  0  THEN  READ  IN  CURRENTS, IF  WK  OT  100  THEN  ZERO  OUT  CURRENTS 
3570  IMPLICIT  DOUBLE  PRECISION  !A-H.  0-Z ) 

3580  PARAMETER  NX-5. NY-6.  N3- (NX+l > *!NY*1 >. NX t-NX  +  l 

3590  Si.  NY1-NY+1.  NXM1-NX-1 .  NYM1  -NY-I ,  N4-NX*NY 


3600  DIMENSION  GAIRX ( NX,  NV > .  GAIRY (NX. NY 1 , GWATXfNX, NY),  CWATY(NX. NY) 

3610  *  FIRST  READ  IN  WIND 
3620  DO  101  1  /  NY 

3630  DO  101  1=1, NX 

3640  GAIRX  (  I .  J)  -*3  OD+OO+O  1  0+00*(  I+J  ) 

3650  GAIRYt  I,  J>=»-4  0D+00*0  2D+00* ( 1+ J > 

3660  GUATX ( I ,  J>«0  01D+00-0.  001D+00*( I-J> 

3670  GWATY(I,J)=0  05D+00+0  001 D+OOM  I  +J  > 

3680  101  CONTINUE 

3790  RETURN 

3800  END 

3810  SUBROUTINE  GROATE ( H, HO, FH, FO, AREA ) 

3820  *  SUBROUTINE  TO  OBTAIN  THICK  AND  THIN  GROWTH  RATES 
3830  IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z> 

3840  PARAMETER  NX=S, NY=6, N3= < NX ♦ 1 > * ( NY+ 1 > , NX  1 -NX  + 1 

3850  fc, NY 1 =NY+ 1 ,  NXM1 =NX- 1 ,  NYM1 =NY- 1 

3860  DIMENSION  HINX1,  NY1 ),  FOtNXl,  NY 1 ) , FH < NX  1 , NY  1 ) 

3870  &, AREA ( NX  1 ,  NY 1,3) 

3830  COMMON/GROWTH/FGROWI21 > 

3890  ONE- 1  O 

3900  ZERO-O  0 

3910  S2-H0+1  0 

3920  M2=IDINT(S2> 

3930  DO  101  J=1,NY1 

3940  DO  101  1=1, NX1 

3950  S1=H( I, U)/0. 5* 1  0 

3960  S1-DMIN1  (SI,  2  0D+01  ) 

3970  Ml  * IDINT  <  S 1  ) 

3980  F0( I.  J ) * ( FGROW ( M2+ 1 >-FGR0W<M2> ) * ( S2-DFL0AT ( M2 ) > +FCROW ( M2 i 

3990  *  THIN  ICE  IS  ALWAYS  TAKEN  TO  BE  ZERO  THICKNESS 
4000  201  F0( I. J)=FGROW( 1 ) 

4010  202  CONTINUE 

4020  FH  (  I,  J  >  =  <  F  GROW  ( M 1  +■  1  >-FGROW(Ml  >  )  *  (  SI  "DFLOAT  ( Ml  )  >+FGROW<Ml  ) 

4030  101  CONTINUE 

4040  DO  102  J= 1 , NY 1 

4050  DO  102  1  =  1, NX  1 

4060  FH( 1,  J)=FH( I,  J ) # AREA ( I .  J.  2 ) +F0 C I .  J > » ( 1  O-AREAI I,  U, ?>  > 

4070  102  CONTINUE 

4080  RETURN 

4090  END 

4100  SUBROUTINE  GROWTH (HEFF.  AREA, FH.  HO, A22 > 

4110  *  SUBROUTINE  TO  CALCULATE  CHANGE  OF  PARAMETERS  DUE  TO  GROWTH 
4120  IMPLICIT  DOUBLE  PRECISION  (A-H, 0-Z> 

4130  PARAMETER  NX-5, NY=6. N3- < NX* 1 ) * ( NY+1 ) . NX  1 -NX ♦ 1 

4140  NY1-NY+1,  NXM1  — NX- 1 .  NYMI=NY-l 

4150  *  THIS  SUBROUTINE  MUST  BE  CALLED  AFTER  ADVECTION 

4160  *  FH  IS  ADDITION  MELTING  THAT  MUST  BE  DONE  TO  BALANCE  MIXED  LAV£R 
4170  DIMENSION  HEFF < NX  1 , NY  1 , 3 > , AREA (NX  1 , NY  1 , 3 > , GHEFF( NX  1 , NY  1 ! 

4180  «. GAREA<  NX  1 , NY1 ) ,  H< NX  1 ,  NY  1 > ,  FH ( NX  1 ,  NY 1 ) 

4190  COMMON/RATE 'FHEFF (NX  1 .  NY  1 >  ,  F0( NX  1 ,  NY l > 

4200  COMMON/STEP/DELTAT, DEL  TAX, DELTAY 

4210  COMMON/ ARRAY/HEFFM( NX  1 . NY 1 ) , UVM(NX, NY) 

4220  COMMON/OUTFLO/QUT ( NX  1  ■  NY 1 > 

4230  ZERO-O  0 

4240  ONE- 10 

4250  *  FIRST  SOLVE  FOR  THICKNESS  OF  THICK  ICE 
4260  DO  101  J— 1 , NY 1 

4270  DO  101  I  —  1 ,  NX  1 

4280  AREA  < I , J.  2 ) -DMAX 1 ( A22,  AREA ( I , J. 2 ) ) 

4290  H( I, J) -HEFF ( I ,  J,  2 1 /AREA ( I ,  V,  2) 

4300  101  CONTINUE 

4310  *  GROWTH  SUBROUTINE  CALCULATES  TOTAL  GROWTH  TENDENCIES  INCLUDING  SNOWFAL 
4320  CALL  OROATE ( H, HO, FHEFF , FO, AREA  > 

4330  1  FORMAT (IX. 11C11.  5) 

4340  #  NOW  CALCULATE  CORRECTED  QROWTH 
4350  DO  201  U— 1 ,  NY  1 
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00  201  1  =  1  >  NX  1 

GHEFF ( I ,  J ) — DELTAT*FHEFK  < 1 .  J > 

GAREA ( I .  J ) -DELTAT*F0 ( I . J i 
201  CONTINUE 

DO  102  J— 1 .  NV 1 
DO  102  1=1, NX1 

GHEFF <  I ,  U>  —  1  0*DMIN1  <HEF'F(  I.  J,  1 ) .  0HEFF (  ]  ,  J >  ) 

GAREA < I, J1-DMAX1 < ZERO. GAREA< 1. J) ) 

102  CONTINUE 

DO  103  J»1 . NY 1 
DO  103  J  "  1  >  NX  1 

FH( X, U)=DMIN1 (ZERO,  GHEFF ( I.  J) ) 

103  CONTINUE 

*  NON  DO  AREA  CHANGE 

DO  104  J=1 . NY1 
DO  104  1*1 , NX  1 

GAREA < I , J>«2  0*( 1.  O-AREA ( 1.  J,  2>  >*GAREA< I,  J) /HO 
«<♦.  5*FH<  I.  J)*AREA<  I,  J,  2  >  /  <  HEFF  <  I ,  J,  1  >  +  00001  > 

104  CONTINUE 

*  NOW  CORRECT  AREA  AND  HEFF<1. J,  11 

DO  105  J= 1 , NY 1 
DO  105  1=1, NX1 

AREA ( I , J,  1 ) “AREA ( I.  J,  1 ) +GAREA ( I ,  J) 

HEFF < I, J, 1 >=HEFF< I, J, 1 ) +GHEFF ( 1 . J)*OUT( I, J) 

105  CONTINUE 

*  NOW  ZERO  OUTSIDE  POINTS 

DO  109  J=1 , NY1 
DO  109  1=1, NX1 

AREA ( I ,  J,  1 ) =AREA < I , J,  1 ) *HEFFM (  I.  J > 

HEFF  < I . J,  1 >=HEFF< I, J.  1 >*HEFFM( I. J) 

109  CONTINUE 

*  NOW  SET  AREA< I,  J,  1 )=0  WHERE  NO  ICE  IS 

DO  110  J= 1 , NY 1 
DO  110  1=1, NX1 

AREA  <  I ,  Ji  1  >=DITIN1  (AREA( I, J.  1  > ,  HEFF  (  I ,  J,  1  )/  0001  ) 

110  CONTINUE 

*  NOW  TRUNCATE  AREA 

DO  106  U=  1  >  NY  1 
DO  106  1=1, NX1 

AREA<I,  J,  1  >=DMIN1  (ONE.  AREAU,  U,  1>> 

106  CONTINUE 

DO  107  J-1.NY1 
DO  107  1=1, NX1 

AREA ( I , J,  1 >  =DMAX 1 ( A22.  AREA ( I , U,  1  >  > 

107  CONTINUE 

*  NOW  CALCULATE  ADDITIONAL  ICE  TO  BE  MELTED  FOR  MIXED  LAYER  BALANCE 

DO  108  J*l» NY1 
DO  108  1=1. NX1 

FH ( I . U ) =CHEFF ( I,  J) -DELTAT*FHEFF < I ,  J ) 

*  NOW  STORE  COMMON  GROWTH  RATE 

FHEFF ( I, J) “GHEFF ( I ,  J) 

108  CONTINUE 
RETURN 
END 

*  MAIN  DRIVING  PROGRAM  FOR  VISCOUS  PLASTIC  SEA  ICE  MODEL 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

*  SPACER 

PARAMETER  NX=5.  NY=6,  N3«(NX-H  )* <NY+1 ) ,  NX  1«NX«-1 
&. NY1-NY+1. NXM1=NX-1,  NYM1=NY-1, N4=NX*NY 
DIMENSION  UICE(NX.  NY,  3>.  VICE(NX, NY,  3) ,  ETA< NX  1 ,  NY1 ) 

ZETA(NX  1 .  NY1  ) ,  DRAOSCNX,  NY).  DRAOAINX,  NY),  OAIRXINX,  NY  > 
it,  0AIRY1NX.  NY ) .  GWATX  (NX,  NY),  QWATY(NX.  NY) ,  HEFF  (NX  1,  NY1 . 3) 
tt.  AMASS  (NX,  NY).  FORCEX(NX,  NY) ,  FORCEY(NX,  NY),  AREA(NX1 ,  NY1 , 3) 
to  HCORR  <  NX 1 ,  NY 1 ) ,  HEFF 1 ( N3 ) 

So  UICEC (NX,  NY),  VICEC(NX,  NY) 

to  0UT1 (NX1 ,  NY1 ) ,  0UT2(NX 1 ,  NY1 >,  UERR(NX. NY), VERR(NX,  NY > . UERRV(N4 > 
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8,, VERRV(N4> 

So  FGRW( 363> 21  > ,  HD IFF < NX  1 .  NYl > ,  FHE1 ( N3 > .  OR  ( N3 ) .  AREA 1 ( N3 ) 
COMMON/RATE/FHEFF (NX  1 ,  NYl ). F0<NX1. NY1 ) 

COMMON/STEP /DELTAT,  DEL TAX,  DELTAV 
C0MM0N/ARRAY/HEFFM(NX1.  NYl), UVM(NX.  NY) 

COMMON/QUTFLO/QUT (NX1, NY 1 > 

CQMMGN/GR0WTH/FGR0W(2l ) 

COMMON/STRESS/STRESS < NX  1,  NY1 , 3) 

COMMQN/DIFF/HDF (NX1, NYl > 

EQU I VALENC  E ( HD I FF , GR  > 

EQU I VALENC  E  <  FHEFF , FHE1  • 

EQU I VALENC E(FHEFF, FHE 1 i 
EQUIVALENCE (AREA. AREA1 > 

EQUIVALENCE  (HEFF, HEFF1 ) 

EQUIVALENCE  ( UERR,  UERRV) .  (VERR.VERRV) 

*  NOW  DECIDE  ON  BASIC  PARAMETERS 

DELTT-86400.  0 
DELTAT*86400.  0 
DELTAX=1  25D+05 
DELTAY*1  25D+05 
ERROR*  00000 1D+00 
ERR0R=5  0D+00*ERR0R 
FHSUM 1 =0  0 
GRSUM1*0  0 
AR SUM 1*0  0 

DIFF1*  002D+00*DEL TAX 
D IFF 1*2  0*DIFF1 
A22*0  15D+00 
H0«0. 5 

*  DOUBLE  HO  BECAUSE  OF  MOD  IN  GROWTH 

H0*2. 0*H0 

LAD=2 

KSTOP*l 

*  NOW  INITIALIZE  COUNTER 

I COUNT *0 
T0UT=0  0 
HESUM*0  O 
UVSUM=0  0 
HHESUM»0  0 

*  NOW  DEFINE  BOUNDARIES 

CALL  BNDRY 

*  WR I TE ( 20 )  ( (OUT( I, J), 1*1, NX1 ). J»l. NYl ) 

DO  122  J*l, NYl 

DO  122  1*1. NX1 

OUT 1 ( I . J ) *  HEFFM ( I, J) -OUT ( I ,  J ) 

HESUM-HESUM+OUT ( I,  J) 

HHESUM*HHESUM*HEFFM( I.  J) 

122  CONTINUE 

DO  125  J«l. NY 
DO  125  1*1,  NX 
UVSUM*UVSUM*UVM( I,  J) 

123  CONTINUE 

PRINT  3.  ( ( HEFFM ( 1, J).  I -1 ,  NX >  ,  J*1 , NY > 

PRINT  4 

PRINT  3.  ( ( UVM ( I ,  J ) ,  I-t, NX), U-l,  NY) 

PRINT  4 

PRINT  3.  ( ( OUT ( I ,  J) ,  1*1 . NX  > ,  U*1 ,  NY  > 

PRINT  4 

PRINT  6. HESUM,  UVSUM 
PRINT  6. HHESUM, UVSUM 

6  FORMAT  (IX.  'HESUM  IS 01  3.  6,  'UVSUM  IS'. 013.7) 

*  NOW  READ  IN  OEOSTROPHIC  WIND  AND  CURRENTS 

*  IF  WK  LT  0  THEN  NEW  CURRENTS  READ  IN,  IF  WK  OT  100  THEN  CUR  ZEROED 

*  31  IS  WIND  FILE 

*  30  IS  CURRENT  FILE 

WK— 1,  OD+OO 
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5580  CALL  CF.O(GAIRX,  OAIPY,  GWATX.  GWATY,  WK  > 

5590  PRINT  3.  ((GAIRXtl.  J),  1-1.  NX>,  J=1,NY> 

5600  PRINT  4 

5610  PRINT  3.  <  (GAIRY<  I,  J),  I  =«  1 ,  NX),  J»l,  NY) 

5620  PRINT  4 

5630  PRINT  3,  ( <GWATX( I,  J).  1  =  1. NX >, J=l. NY) 

5640  PRINT  4 

5650  *  NOW  SET  UP  GROWTH  RATES 

5660  CALL  FGROWP 

5670  *  FILE  32  CONTAINS  GROWTH  RATES 

5680  *•  READ  (32)  f  Gf, W 

5690  *  NOW  INITIALIZE  SYSTEM 

5700  *  FIRST  GUESS  AT  INITIAL  HEFF  AND  AREA 

5710  DO  101  J= 1 , NY 1 

5720  DO  101  1=1, NX1 

5725  HDIFF  < I ,  J ) =0. 0 

5730  HEFF ( I , J,  3 ) =0  0 

5740  HEFF ( I . U,  2)=0  0 

5750  AREA ( 1 1 J.  2 )  =  1  0 

5760  AREA< I- J,  3>=1.  0 

5770  HEFF ( I,  J.  1  )»(3. OD+OOl/O.  91D+00 

5780  HEFF ( I , J,  1 ) =HEFF < 1 , J,  1 ) *0U7  < I ,  J) 

5790  AREA( I. J. 1 >=1. 0 

5800  101  CONTINUE 

5810  CALL  XSUM(HEFF1 , THEFF) 

5820  *  NOW  ZERO  OUT  VELOCITY 
5830  DO  102  J=  1 1  NY 

5840  DO  102  1=1, NX 

5850  UICE< I, J, 3)=0.  0 

5860  UICE ( I , J,  2 ) =0  0 

5870  UICE < I, J,  1 )=0.  0 

5880  VICE < I , J, 3 ) =0.  0 

5890  VICE ( I , J,  2  >=0.  0 

5900  VICE ( I. J,  1 )=0.  0 

5910  UICECd,  J)=0.  O 

5920  VICEC ( I , J ) =0  O 

5930  102  CONTINUE 

5940  *  NOW  GET  FIRST  VALUE  OF  U  AND  V 
5950  THETA=1  O 

5960  CALL  FORM(UICE.  VICE. ETA.  ZETA. DRAGS. DRAGA, CAIRX, CAIRY, CWATX, GWATY 

5970  FORCEX,  FORCEY,  HEFF.  AMASS,  AREA) 

5980  PRINT  3,  < (FORCEX  < I,  J) ,  I  =  1 .  NX  > ,  J-l , NY ) 

5990  PRINT  4 

6000  PRINT  3.  < (FORCEY( I, J),  1  =  1,  NX), J-l. NY) 

6010  *  NOW  SET  AMASS-0  AND  DEFINE  INITIAL  VISCOSITY 
6020  DO  103  U=1,NY 

6030  DO  103  1=1, NX 

6040  AMASS ( I.  J)=0.  0 

6050  ZETA(  I,  J)=HEFF( I,  U.  1  )*(  1  OD+U  > 

6060  ETA( I,  J)=ZETA( I, J)/4  0 

6070  103  CONTINUE 

6080  CALL  RELAX (UICE. VICE, ETA,  ZETA,  DRAGS, DRAGA, AMASS, FORCEX,  FORCEY 

6090  «<,  ERROR,  THETA,  UICEC,  VICEC  ) 

6100  PRINT  4 

6110  PRINT  4 

6120  PRINT  3,  ( <UICE( I, J.  1 ).  1  =  1. NX), J-l. NY) 

6130  PRINT  4 

6140  PRINT  3,  ( (VICE( I. J.  1 ).  1-1, NX), J-l. NY) 

6150  PRINT  4 

6160  PRINT  4 

6170  *  NOW  START  STANDARD  PREDICTOR  CORRECTOR  PROCEDURE 

6180  100  CONTINUE 

6190  *  NOW  PUT  GROWTH  RATES  IN 

6200  KORO-MOD( I COUNT, 365)+l 

6210  «  DO  141  KCO-1, 21 

6220  •  FOROW( KGO ) >FQRW( KORO,  KCG ) / <  8  64D+06 ) 
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*141  CONTINUE 

*  NOW  SET  UP  REASONABLE  OUTFLOW  THICKNESSES 

CALL  MEAN ( HEFF ,  0UT2  > 

DO  123  J-2, NV 
DO  123  1-2,  NX 

HEFF  <  I .  J.  1)=HEFFCI. J.  I  >  +QUT1  <  I ,  J)*0UT2U ,  J> 

123  CONTINUE 

CALL  MEAN ( AREA,  0UT2) 

DO  124  J-2, NY 
DO  124  I =2, NX 

AREA! I, J,  1 )=AREA( I, J,  1 )>0UT1 ( I,  J)*0UT2C I,  J) 

124  CONTINUE 

CALL  X  SUM  <  HEFF 1 , THEFF 1 ) 

THEFF  t -THEFF 1 -THEFF 

*  FIRST  DO  PREDICTOR 

DO  121  J=1,NY 
DO  121  1=1, NX 
UICE< 1 , 3, 3>=UICE < I,  J,  1) 

VICE! I, J,  3>=VICE< I,  J,  1) 

UICECCI,  J)  =UIC£ ( I .  J,  1) 

VICECCI. J)=VICECI.  J.  1) 

121  CONTINUE 

*  DO  121  J=1 , NY 

*  DO  121  1=1, NX 

*  UICECU,  J)«UICE(I,  J,  1> 

*  VICEC (I,  J) -VICEC I,  J,  1 ) 

*121  CONTINUE 

THETA- 1 .  0 

*  DELTAT-DELTT/2  0 
DELTAT-DELTT/2  0 

CALL  FORMCUICE.  VICE,  ETA, ZETA.  DRAOS. DRAOA. GAIRX,  DAIRY, CWATX. GWATY 
&.  FORCEX,  FORCEY,  HEFF.  AMASS.  AREA) 

CALL  RELAX  CUICE, VICE. ETA, ZETA, DRAGS. DRAGA, AMASS.  FORCEX.  FORCEY 
ERROR.  THETA,  UICEC,  VICEC) 

*  NOW  DO  REGULAR  TIME  STEP 

*  NOW  DOING  BACKWARDS  TIME  STEP 

THETA- 1  O 
DELTAT-DELTT 

CALL  FORMtUICE,  VICE,  ETA,  ZETA. DRAGS, DRAGA,  GAIRX.  DAIRY,  GWATX,  GWATY 
&,  FORCEX,  FORCEY,  HEFF,  AMASS,  AREA) 

*  NOW  SET  U< 1 >=U<2>  AND  SAME  FOR  V 

DO  111  J-l.NY 
DO  111  1  =  1 ,  NX 
UICEU.  J.  3)=UICE<  I,  J.  1  ) 

VICEU.  J.  3)-VICE<  I,  J,  1  ) 

UICEC ( I.  J) -UICEC I, J, 1) 

VICEC ( I. J> -VICE (1,3, 1) 

*  UICEC <  I.  J ) -UICEC  I.  J,  1 ) 

*  VICEC < I. J) -VICEC I. J, 1 ) 

UICEC I. J. 1 ) -UICEC I, J.  2) 

VICEC I. J,  1 > -VICEC 1,  3, 2) 

111  CONTINUE 

CALL  RELAX  CUICE,  VICE.  ETA,  ZETA.  DRAOS,  DRAGA,  AMASS,  FORCEX.  FORCEY 
«t.  ERROR.  THETA,  UICEC,  VICEC  ) 

ICOUNT-ICOUNT+1 

*  NOW  DO  ADVECTION 

CALL  ADVECT CUICE,  VICE,  HEFF,  DIFF1,  LAD) 

CALL  ADVECT  CUICE,  VICE.  AREA,  DIFF1,  LAD) 

«  NOW  DO  GROWTH 

CALL  QROWTHCHEFF, AREA,  HCORR,  HO,  A22) 

*  MUST  CALL  GROWTH  ONLY  AFTER  CALL I NO  ADVECTION 

*  NOW  CORRECT  OUTFLOW  PTS  AND  GET  OUTFLOW  ICE 

CALL  XSUM C  HEFF 1 ,  THEFF ) 

DO  105  J-t.NYl 
DO  105  I-1.NX1 

HEFF ( I . J.  1 >-HEFF ( I.  J.  1 > *OUT ( I ,  J > 
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6880  AREA  < I ,  J,  1  )  “AREA  <  I ,  J,  1  )  #CJUT  <  I .  J  > 

6885  FHEFF! I,  J ) =FHEFF ! I , J ) *0UT ! I ,  J) 

6886  HDIFF< I. U)  =  !  1  0-AREA < I ,  J,  2) )*F0< I,  J)*0UT < I,  J)*DELTT 

6890  105  CONTINUE 

6900  CALL  XSUM  !  HEFF 1 ,  THEFF2  > 

6902  CALL  XSUM ( AREA 1 ,  ARSUM ) 

6903  CALL  XSUM  <  FHE l ,  FHSUM ) 

6904  CALL  XSUM (GR<  GRSUM) 

6905  OR  SUM 1  =GR  SUM 1 + GR  SUM 

6906  FHSUM 1 =FHSUM 1+FHSUM 

6907  ARSUM1 “ARSUM1 +ARSUM 

6910  TOUT 1 “THEFF-THEFF2-THEFF 1 

6920  THEFF=THEFF2 

6930  TOUT=TOUT+TOUT 1 

6935  ICOU=ICOUNT 

6940  *  NOW  PRINT  OUT  EVERY  KTH  PT 

6945  *  IF! ICOUNT.  LT.  366)  00  TO  97 

6946  *  I COU= I COUNT-365 

6950  KTH=MOD< ICOUNT,  8) 

6955  KTH=MOD( ICOUNT, 20) 

6960  IF1KTH.  EQ  1)  GO  TO  95 

6970  GO  TO  96 

6980  95  CONTINUE 

6990  4  FORMAT <//) 

7000  1  FORMAT! IX,  'TIME  STEP  AND  TOTAL  THICKNESS  ARE',2G20.  12) 

7010  PRINT  1 ,  ICOUNT, THEFF1 

7020  PRINT  1,  ICOUNT, THEFF 

7030  PRINT  4 

7040  PRINT  2,  T0UT1,  TOUT 

7050  PRINT  4 

7060  2  FORMATUX,  'OUTFLOW  AND  NET  OUTFLOW  ARE'.2C14  7) 

7070  PRINT  3,  < (UICE< I, J,  1 ),  1  =  1 , NX ) ,  J*1 ,  NY ) 

7080  PRINT  4 

7090  PRINT  3,  < <VICE( I,  J,  1 ),  1*1,  NX),  J»l,  NY) 

7100  PRINT  4 

7110  PRINT  3,  ( < HEFF < I , J,  1),  1  =  1 ,  NX ) .  J=1 ,  NY ) 

7120  PRINT  4 

7130  PRINT  3,  (1AREAU,  J.  1),  1*1 ,  NX  >  ,  U-l ,  NY  » 

7140  PRINT  4 

7150  PRINT  3,  !!ZETA!I,  J),  1*1, NX),  J-l.NY) 

7160  PRINT  4 

7170  PRINT  3,  < !FHEFF! I, J),  1*1, NX),  J*l. NY) 

7180  PRINT  4 

7190  PRINT  3,  < (STRESS! I, U.  1 ),  1-1,  NX),  U«l,  NY) 

7200  PRINT  4 

7210  PRINT  3,  <FGROW( I ),  1*1, 21 ) 

7220  PRINT  4 

7230  3  FORMAT! IX, 5G20  14) 

7240  *  NOW  PRINT  OUT  5QUARE  VELOCITY  AND  VELOCITY  DIFERENCE  AND  MAX  CHANCE 

7250  SQ1=0  0 

7260  SQ=0. 0 

7270  DO  130  U=  1 ,  NY 

7280  DO  130  1=1. NX 

7290  SQ=SQ+UICE ( I ,  J,  1 )**2+VIC£< I,  J,  1 >**2 

7300  UERR < I , J)=UICE( I.  J,  1 )-UICE( I, J. 2) 

7310  VERR ( I , U)=VICE< I, J,  1 >-VICE< I.  J,  2> 

7320  SQ 1 =SG 1 +UERR ( I , J) *#Z+VERR < I , U>**2 

7330  130  CONTINUE 

7340  *  ISU=MAXMAC (UERRV ) 

7350  #  ISV-MAXMAG(VERRV) 

7355*  SMU=DABS< UERRV! ISU+1 ) ) 

7356  *  SMV=DABS  <  VERR V ( I SV+ 1 ) ) 

7360  CALL  XMAXM ! UERRV. SMU) 

7370  CALL  XMAXM ! VERRV, SMV) 

7380  SM-DMAX 1 ! SMU.  SMV ) 

7390  PRINT  5. SQ, SQ1.SM 


7400  5  FORMAT (IX.  'SQUARE  VELOCITY,  SQ.  VELOCITY  DIFFERENCE.  MAX  CHANCE 

7410  S./1X.3G20  12) 

7420  PRINT  4 

7430  96  CONTINUE 

7440  *  WRITE! 20 )  < < UICE < I , J,  1 ) ,  1*1 , NX ) . U-l , NY ) 

7450  *  WRITE! 20)  ( ( VICE ( I ,  J,  1 ) ,  1  =  1,  NX > , J-l , NY ) 

7460  *  WRITE  (20)  < ( HEFF <I,J,  1),  I  =  1,NX1),J=1,NY1) 

7470  *  WRITE  (20)  ( ( AREA ( I , J,  1 ) ,  1-1 , NX  1 ) , J»1 , NY1 ) 

7480  *  WRITE  (20)  ( (FHEFF ( I,  J,  1), 1  =  1 ,  NX1 > ,  J=1 , NY  1 > 

7490*  WRITE  (20)  < ( ZETA< I, J) , 1=1, NX1 ), J-l, NY1 ) 

7500  *  WRITE  (20)  ( ( STRESS < I, J,  1 ) ,  1*1 , NX  1 ) , J=1 , NY1 ) 

7505  *  WR I TE ( 23 )  TOUT 1 , FHSUM,  GR SUM, ARSUM, THEFF 

7510  97  CONTINUE 

7511  PRINT  1,  ICOU, THEFF 

7512  PRINT  2, TOUT 1 ,  TOUT 

7513  PRINT  8, FHSUM,  FHSUM 1 

7514  PRINT  9, GRSUM, GRSUM1, ARSUM,  ARSUM1 

7515  9  FORMAT (IX,  'OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC',4G15  7) 

7516  8  FORMAT (IX,  'GROWTH  AND  NET  GROWTH  ARE',2G14.  7) 

7520  *  NOW  CHECK  COUNTER  AND  DECIDE  IF  NEW  WINDS  NEEDED 

7530  KTH1=M0D< ICOUNT,  365) 

7540  IF(KTH1EQ  0)  GO  TO  211 

7550  GO  TO  212 

7560  211  CONTINUE 

7570  *  REWIND  31 

7580  PRINT  7 

7590  7  FORMAT  (IX, 'WIND  FILE  REWOUND') 

7600  212  CONTINUE 

7610  *  READ (31 )  GAIRX, GAIRY 

7620  *  NOW  DECIDE  IF  DONE 

7630  IF ( ICOUNT.  GT  50)  GO  TO  200 

7635  IF( ICOUNT  GT.  20)  GO  TO  200 

7640  GO  TO  lOO 

7650  200  CONTINUE 

7660  KSTOP— 2 

7670  IF  (KSTOP.  EQ  1)  GO  TO  202 

7680  GO  TO  201 

7690  202  CONTINUE 

7700  DO  204  J=1,NY 

7710  DO  204  1=1, NX 

7720  GAIRX ( I,  J)=0.  0 

7730  GAIRY ( I,  J)=0.  O 

7740  204  CONTINUE 

7750  IC0UNT=0 

7760  KST0P=2 

7770  GO  TO  100 

7780  201  CONTINUE 

7790  *  WR I TE ( 2 1 )  UICE, VICE,  HEFF, AREA 

7800  STOP 

7810  END 

7820  SUBROUTINE  ME AN (HEFF, HME AN) 

7830  *  SUBROUTINE  FINDS  MEAN  HEFF  AT  OUTFLOW  PTS  OF  VALUES  AROUND 
7840  IMPLICIT  DOUBLE  PRECISION  (A-H, 0-Z) 

7850  PARAMETER  NX  =  5, NY=6. N3= (NX+l >* (NY+l ) >  NX l=NX+l 

7860  NY1=NY+1,  NXM1=NX-1.  NYM1 *NY- 1 ,  N4=NX*NY 

7870  DIMENSION  HEFF (NX  1 , NY1 , 3 ) , HMEAN (NXI, NY1 > 

7880  COMMON/OUTFLO/OUT ( NX  1 , NY1 > 

7890  DO  101  U=2, NY 

7900  DO  101  1=2, NX 

7910  HMEAN ( I >  J ) = ( HEFF ( 1  +  1, J,  1 >*OUT< 1  +  1.  J)*HEFF( 1*1 , J*1 . 1 >*OUT(I  +  l, J+l ) 

7920  $,  +HEFF( 1*1, J-l, 1 >*OUT( 1  +  1, J-l >+HEFF( I, J+l, 1 >*OUT( I, J+l ) 

7930  &  +HEFF ( I , J-l ,  1 ) *OUT ( I ,  J-l ) +HEFF ( 1-1 ,  J,  1 ) *OUT ( 1-1,  J) 

7940  S<  +HEFF< 1-1, J+l,  1 >*OUT( I-i,  J+l >+HEFF( 1-1, J-l.  1 >*OUT( 1-1,  J-l ) 

7950  t< ) / ( OUT ( 1  +  1,  J)+OUT( 1+1,  J+l >+OUT( 1  +  1, J-l >+OUT( I, J* 1 ) +OUT ( I , J-l ) 

7960  *  +OUT ( 1-1 , J)*OUT ( 1-1 ,  J+l )+OUT ( 1-1 , J-l >  +  .  00001D+00) 

7970  101  CONTINUE 


7980 
7990 
8000 
8010 
8020 
8030 
8040 
8050 
8060 
8070 
8080 
8090 
8100 
81  10 
8120 
8130 
8140 
8150 
8160 
8170 
8180 
8190 
8200 
8210 
8220 
8230 
8240 
8250 
8260 
8270 
8280 
8290 
8300 
8310 
8320 
8330 
8340 
8350 
8360 
8370 
8380 
8390 
8400 
8410 
8420 
8430 
8440 
8450 
8460 
9000 
9010 
9020 
9030 
9040 
9050 
9060 
9070 
9080 
9090 
9100 
9110 
9120 
9130 
9140 
9150 
9160 


RETURN 

END 

SUBROUTINE  PLASTIUICE, VICE, PRESS, ETA, ZEYA, ECCEN, ZMAX, ZMIN) 

►  SUBROUTINE  CALCULATES  STRAIN  RATES  AND  VISCOUS  PARAMETERS 
IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z  ) 

PARAMETER  NX=5, NY=6,  N3= ( NX+1 > * ( NY+1 > ,  NX1=NX+1 
!<,  NY  1  =NY+ 1 ,  NXM1  =NX  — 1 ,  NYM 1  —NY—  1 ,  N4=NX*NY 
DIMENSION  U  ICE  (NX,  NY,  3) ,  VICE  (NX,  NY,  3) ,  PRESS  ( NX  1.  NYU 
ETA< NX  1 ,  NY  1  ) ,  ZET A ( NX  1 ,  NY!  ) ,  El  1  (NX  1  ,  NY1  ) ,  E22 ( NX  1 ,  NY1  ) 

So  E12 ( NX  1 , NY1 ) ,  ZMAX ( NX  1, NY 1  ) ,  ZMI N ( NX  1 ,  NY 1 ) 

COMMON/STEP/DELTAT,  DELTAX,  DEL.TAY 
COMMON/STRESS/STRESS (NX  1 ,  NY  1 , 3) 

COMMON/ARRAY/HEFEM( NX  1 ,  NY  1 ) ,  UVM(NX,  NY) 

ECM2=1  0/ <  ECCEN**2 ) 

GMIN=1  0D-20 

*  NOW  EVALUATE  STRAIN  RATES 

DO  101  U-2, NY 
DO  101  1=2, NX 

El  1 ( 1, J)  =  (0  5/DELTAX )*(U1CE( 1, J,  1 >+U:CE( I,  J-l.  1 ) 

Z,  -UICE<  1-1,  J,  1  )-UICE(  1-1,  J-l,  1  )  ) 

E22< I, J>  =  (0  5/ DELTA Y ) * ( V I CE ( I , J, 1 ) +V I CE ( I - 1 , J. 1 > 

8<  -VICE <  I,  J-l,  1  > -VICE ( 1-1,  J-l,  1  >  > 

E12( I, J  >  =•  ( 0 .  25 /CELT AY ) * ( UI CE ( I, J,  t >+UlCE ( 1-1 , J,  1 ) 

&  — UlCE (  I ,  J-l,  1  )-UICE(  l- 1,  J-l,  1  >  > 

?<  +(0  25/DELTAX  )  *  <  VICE  (  I ,  J,  1  )  WICE  (  I ,  J-l ,  1) 

&  -VICE ( I  — 1,  J,  1  ) -VICE ( I -1,  J-l,  1 >  > 

*  NOW  EVALUATE  VISCOSITIES 

DELT=(E1 1 ( I,  J ) **2+E22 ( I ,  J>**2>*< 1  0+ECM21+4.  0*£CM2*El2( I , J>**2 

5,  +2.  0*E1 1  (  I,  J)*E22<  I,  J>*(  1.  0-ECM2) 

DELT 1 =DSQR  T ( DELT ) 

DELT1=DMAX1 ( GMIN,  DELT 1 ) 

ZETA ( I ,  J)=0  5*PRESS(I, J)/DELT1 

101  CONTINUE 

*  NOW  PUT  MIN  AND  MAX  VISCOSITIES  IN 

DO  102  J=l, NY 1 
DO  102  1=1, NX1 

ZETACI,  J)=DMIN1!ZMAX(I,  J>,  ZETAU,  J>> 

ZETA ( I ,  J ) =DMAX 1 ( ZMIN< I , J) > ZETA ( I , J ) > 

ETA ( I , J)=ECM2*ZETA( I,  J) 

Ell! I,  J)=E11 < I,  J)*HEFFM( I, J> 

E22 ( I , J ) =E22 < I , J  >  *HEFFM  < I , J ) 

E12(I,J)=E12(I,  J)*HEFFM! I.  J) 

SSI  1  =  (  ZETA (  I,  J)-ETA<  I,  J>  )*(El  1(1,  J)-*E22(  I,  J)  ) -PRESS!  I,  J)*0.  5 
STRESS! I, J,  1 )  =  <2.  0*ETA( I, J)*E1 1 ( I,  J)+SS1 1 ) 

STRESS! I, J, 2)=2.  0*ETA! I,  J>*E22( I,  J)+SS1 1 
STRESS! I,  J, 3 >  =2  0*ETA( I,  J)*E12( I,  J) 

102  CONTINUE 
RETURN 
END 

SUBROUTINE  RELAX !OUICE,  OVICE, OETA,  QZETA, ODRAGS, ODRAGA 
S<,  OAM ASS,  ORCEX,  ORCEY,  ORROR,  OTHETA,  OUICEC,  OVICEC) 

IMPLICIT  DOUBLE  PRECISION  !A-H, 0-Z) 

PARAMETER  NX=5, NY=6, N3=(NX+1 >* < NY+1 > ,  NX  1 -NX  +  1 
S<,  NY1-NY+1,  NXM1=NX-1,  NYM1  =NY- 1 ,  N4=NX#NY 
DIMENSION  UICE ! NX, NY, 3), VICE1NX, NY, 3 > , ETA! NX  1 , NY1 ) .  ZETA!NX 1 ,  NY1 ) 

6,  DRAGS! NX,  NY), DRAGA ( NX ,  NY > ,  FORCEX!NX,  NY),  FORCEY ! NX,  NY),  FXETA(4) 
«>,  FXZETA !  4 ) ,  FYETA !  4  ) ,  FYZETA  (  4  ) ,  UERR  ( NX ,  NY ) ,  UERRV  !N4 ) 

&,  VERRINX, NY).  VERRV(N4),  COEF!NX,  NY) . AMASS! NY- NY) 

«<,  FXMINX,  NY),  FYM!NX,  NY) 

So  UICEC ! NX, NY),  VICEC!NX,  NY) 

5,  FXE<4.NX.NY).FYE!4,NX.NY)>FXZ(4,NX,NY>,FYZ(4,NX,NY) 

6,  COEFI!NX,NY) 

&, OUICE!NX, NY, 3), OVICE (NX, NY, 3 ) , OETA (NX  1 >  NY1 ) , OZETA(NX 1 , NY1 ) 

&, ODRAGS (NX, NY), ODRAGA (NX, NY ) , OAMASS (NX, NY),  ORCEXINX,  NY) 

S, ORCEY (NX.  NY),  OUICEC(NX,  NY ) ,  OVICEC ( NX ,  NY ) 

*<,  HEFFM(NX  1 ,  NY1  ) ,  UVM!NX,  NY) 
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9170 

9180 

9190 

9200 

9210 

9220 

9230 

9240 

9250 

9260 

9270 

9280 

9290 

9300 

9310 

9320 

9330 

9340 

9350 

9360 

9370 

9380 

9390 

9400 

9410 

9420 

9430 

9440 

9450 

9460 

9470 

9480 

9490 

9500 

9510 

9520 

9530 

9540 

9550 

9560 

9570 

9580 

9590 

9600 

9610 

9620 

9630 

9640 

9650 

9660 

9670 

9680 

9690 

9700 

9710 

9720 

9730 

9740 

9750 

9760 

9770 

9780 

9790 

9800 

9810 

9820 


COMMON/D 1NV/DELIN2 

COMMON/STEP/OTAT,  OTAX,  OTAV 

COMMON /ARRAY/OHEFFM< NX  1.  NY1 ) ,  OUVMCNX,  NY) 

EQUIVALENCE  (UERR.  UERRV) , (VERR. VERRV) 

IC0UNT=0 

WFA=1.  5  _ 

*  NOW  REPLACE  DOUBLE  PRECISION  INPUT  BY  SINGLE  PRECISION  VARIABLES 

DO  601  J=l, NY 
DO  601  1  =  1.  NX 
UVM< 1 ,  J ) =OUVM ( I < J) 

DRAGS < I > J)=ODRAGS< I,  U> 

DR AG A < I. J ) “ODRAGA ( I, J> 

AMASS ( I , J)=OAMASS< I,  U> 

FORCEX  < I .  J)=ORCEX< I, U) 

FORCEY  < I . J)=ORCEY( I, J) 

UICEC  < I . J)=OUICEC< I. J) 

VICEC ( I . J)=OVICEC< I. J) 

601  CONTINUE 

DO  603  K=l,  3 

DO  603  J=l. NY 

DO  603  1=1. NX 

UICE ( I > J. K)=0UICE< I, J,  K> 

VICE ( I . J. K >=OVICE< I.  U,  K) 

603  CONTINUE 

DO  602  J= 1 >  NY 1 
DO  602  1=1, NX1 

*  HEFFM(  I,  J)=OHEFFM(I,  J) 

ETA( I,  J)=OETA< I,  J) 

ZETA ( I , J ) =OZETA ( I . J ) 

602  CONTINUE 
ERROR =ORROR 
THET  A=OTHETA 
DELTAX=OTAX 

*  DELT  AY =0T  AY 
DELTAT=OTAT 
DELIN  =1 .  O/DELTAX 
DELIN2=0. 5/ ( DELTAX**2> 

K=1 

*  MUST  UPDATE  HEFF  BEFORE  CALLING  RELAC 

*  FIRST  SET  U ( 2 ) =U < 1 ) 

DO  99  J=l,  NY 
DO  99  1  =  1,  NX 

*  NOW  MAKE  SURE  BDRY  PTS  ARE  EQUAL  TO  ZERO 

UICE  < I, J, 2)=UICE< I, U. 1 > 

VICE ( I , J, 2)=VICE< I. U, 1 ) 

OUICE ( I , J, 2 ) =OUICE  < I , J,  I ) 

OVICE ( I , U,  2)=0VICE< I,  J,  1 > 

OUICE ( I ,  J,  1 )=OUICE( I,  J, 3>*0UVM< I,  J) 

OVICE< I,  J,  1 >=OVICE( I, J, 3>*0UVM< I,  J) 

UICE ( I , J,  1 >=UICE< I, J. 3>*UVM< I.  J> 

VICE ( I , J,  1 )=VICE< I,  J, 3)*UVM( I,  U) 

99  CONTINUE 

*  NOW  SET  UP  COEFFICIENTS  OF  DIAGONAL  COMPONENTS 
501  CONTINUE 

DO  102  U=l, NY 
DO  102  1=1. NX 

COEF< I,  J) “AMASS ( I.  J)/DELTAT+2.  0*THETA*(0.  5*DRAQS< I. J) 

«<+ 2.  0*<  (ETA <  I ,  J) +ETA<  1  +  1  <  J)+ETA(  I.  J+l  >+ETA<  I+l,  J+l  >  > 

&+.  5*(ZETA(  I,  J)+ZETA(I+1, J>+ZETA( I,  J+l  >+ZETA<  1  +  1,  J+l  >  ) 

«<  )  /  <  4 .  0*  ( DELT  AX  **2  >  >  ) 

COEFI < I ,  J ! =1 .  O/COEF ( I,  J) 

102  CONTINUE 

*  NOW  CALCULATE  ALL  FUNCTIONS  OF  PREVIOUS  U  AND  V  VALUES 

TTHETA-2.  0*(1.  O-THETA) 

DO  111  J=2, NYM1 
DO  111  1=2, NXM1 


23 


9830  502  CONTINUE 

9840  CALL  FELLIP (UICE, VICE. ETA. FXETA.  I,  J. 2) 

9850  CALL  FELLIP (UICE,  VICE,  ZETA,  FXZETA,  I , J, 2 ) 

9860  CALL  FELLIP (VICE, UICE. ETA, FYETA, I , J, 2) 

9870  CALL  FELLIP (VICE, UICE, ZETA,  FYZETA,  I , J, 2 > 

9880  FX0=0.  5* (FXETA (I >+FXZETA( 1 ) +FXETA ( 2 )  +f-  XETA( 3) +FXZETA( 4) -FXETA ( 4 ) ) 

9890  FX0=TTHETA*FX0 

9900  FX1=( AMASS ( 1,  J) /DELTAT-TTHETA*0.  5*DRAGS( I, J>  >*UICE( I,  J,  2) 

9910  FX2=TTHETA*0.  5*PR AGA ( I , J)#VICE( I, J, 2) 

9920  FY0=0  5*(FYETA(  1  )+FYETA(  2  )  +FYZETA(2 )  +FYZETA(  3 )  -FV£TA  (  3 )  +FYETA(  4  )  ) 

9930  FY0=FY0*TTHETA 

9940  FY1=  < AMASS ( I , J ) /DELTAT-TTHETA*0  5*DRAGS( I, J> >*VICE( I,  J,  2) 

9950  FY2=-TTHETA*0.  5*DRAGA( I,  J) *UICE ( I , J,  2 ) 

9960  FXC=AMASS ( I , J ) *0.  5*TTHETA* 

9970  St  (UICEC(  I,  J)*(UICE(  1  +  1,  J.  2>-UICE(l-t.  J,  2)  ) 

9980  St  +VICEC  (I,  J)*(UICE(  I,  J+l .  2 ) -UICE  (I,J-l,2>>>/(2.  O+DELTAX  ) 

9990  FXM( I, J )=FX0+FX 1+FX2+FQRCEX ( I , Jl+FXC 

lOOOO  FYC=AMASS< I, J>*0.  5+TTHETA* 

10010  St  (UICEC  (I,  J )* ( VICEl  1  +  1,  J,  2>-VICE(  I  -1 ,  J,  2)  ) 

10020  S<  +VICEC( I, J)*(VICE( I,  J+l, 2)-VICE( I, J-l. 2) ) )/(2.  0*DELTAX ) 

10030  FYM( I, J>=FYO+FYl+FY2+FDRCEY( I,  Jl+FYC 

10040  111  CONTINUE 

10050  *  NOW  SET  U( 3)=U(1) 

10060  100  CONTINUE 

10070  DO  101  J— 1 ,  NY 

10080  DO  101  1=1, NX 

10090  UICE< I , J, 3) =UICE ( I , J, 1 ) 

10100  VICE( I,  J,  3)=VICE( I, J,  1 ) 

10110  101  CONTINUE 

10120  *  NOW  BEGIN  SWEEP 

10130  CALL  FELLD1 (UICE, VICE,  ETA, FXE,  1 ) 

10140  CALL  FELLD1 (UICE, VICE, ZETA, FXZ, 1 ) 

10150  CALL  FELLD1 (VICE,  UICE,  ETA, FYE,  1 > 

10160  CALL  FELLD1 (VICE. UICE, ZETA, FYZ. 1) 

10170  DO  103  J=2, NYM1 

10180  DO  103  1=2,  NXM1 

10190  504  CONTINUE 

10200  K=1 

10210  FXETA(l)  =  FXE< 1,  I,  JJ+DELIN2* 

10220  S<(UICE(  1-1,  J.  R)*(  ETA(I,J+1>+  ETA(I,  J))) 

10230  FXETA(2)=FXE(2, I, J)+DELIN2* 

10240  St  (UICE  (I,  J-l,  KHMETAd.  J)+ETA(I  +  1,  J>  >  ) 

10250  FXETA(3)=FXE(3>  I, J)+0.  5#DELIN2* 

10260  S,  ( VICE  (  I  — 1 ,  J-l,  K)*ETA(  I,  J )  +VICE  (  I ,  J-l ,  K>* 

10270  S<(-ETA( I,  J >+ETA( 1  +  1 , J) )-VICE( 1  +  1.  J-l , K )*ETA ( 1  +  1 ,  J) 

10280  St+VICE(  I  —  1 ,  J,  K)#(-ETA(  I,  J+l  )  *-ETA(  I,  J)  ) 

10290  St-VICE(  1-1.  J+l,  K)#ETA(  I,  J+l  )  > 

10300  FXETA ( 4 >  =FXE ( 4,  I , J ) +DEL IN2*0,  5* 

10310  S<  ( VICE  (  I  — 1 ,  J-l.  K)*ETA(  I,  J)+VICE(  I ,  J-l ,  K>* 

10320  St(-ETA(  1  +  1,  J)+ETA(  I.  J )  ) -VICE ( 1  + 1 ,  J-l,  K>*ETA<  I  + 1 ,  J) 

10330  S<+VICE(  I  — 1 ,  J.  K)*(ETA(  I,  J+l  )-ETA(  I.  J)  ) 

10340  St-VICE(  1-1.  J+l,  K)*ETA(  I,  J+l  >  ) 

10350  FYETA( t )=FYE( 1, I, J1+DELIN2* 

10360  St(VICE(  I  —  1 .  J,  K>*(  ETA(  I ,  J+l  )+  ETA(I,J>>> 

10370  FYETA(2)=FYE(2, I, J1+DELIN2* 

10380  S<  ( VICE ( I,  J-l,  K)*(ETA(  I.  J) +ETA(  1  +  1 .  J )  )  ) 

10390  FYETA(3>=FYE(3, I, J)+0  5+DELIN2* 

10400  St(UICE(  1-1,  J-l,  K)*ETA(  I,  J )  +UICE (  I ,  J-l,  K)* 

10410  St(-ETA( I, J)+ETA( 1  +  1 ,  J ) ) -UICE ( 1  +  1 ,  J-l .  K)»ETA< 1  +  1,  J) 

10420  St+UICE(  I -1 ,  J,  K)*(-ETA(  I,  J+l  >+ETA(  I,  J)  ) 

10430  St-UICE(  1-1,  J+l,  K>*ETA(  I,  J+l  )  ) 

10440  FYETA ( 4 >  =FYE ( 4,  I, J)+DELIN2*0.  5* 

10450  St ( UICE (  1-1,  J-l.  K)*ETA(  I,  J)  +UICE ( I ,  J-l,  K>* 

10460  St ( -ETA(  I  + 1 .  J )  +ETA (  I ,  J)  ) -UICE (  1+ l.  J-l.  K)*ETA(  I  + 1 ,  J) 

10470  St+UICE(  I  — 1 ,  J,  K  )  *  (ETA(  I.  J+l  )-ETA(  I,  J)  ) 

10480  St-UICE ( I- 1 ,  J+l.  X)*ETA(  l,  J+l  )  > 
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10490 
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FXZETAd)-FXZ<l,  I>  J)+DELIN2* 

*c(UICE<  1-1.  J.  K)*(ZETA(I,  J+l  )+ZETA<  1.  J) ) ) 

FXZETA<4)-FXZ<4,  I.  J1+DELIN2+0.  9# 
fc(VICE<  1-1.  J-l.  K)*ZETA<  I,  J)+VICEd.  J-1.  K># 
fc(-ZETA<  1  +  1.  J)+Z£TA<  I.  J>  >-VICE<  1+1.  J-l.  K)#ZETA<  1  +  1.  J> 

S<+VICE(  1-1.  J.  K)#<ZETA<  I.  J+l >-ZETA<I,  J>  ) 

«.-VICE<  1-1.  J+l.  K)*ZETA<  X,  J+l )  > 

FYZETA<2)-FYZ<2,  I. J)+DELXN2» 
fc( VICE(  I.  J-l.  K)«(ZETAd.  J)+ZETA<  1  +  1,  J)  ) ) 

FYZETA<  3)»FYZ (3,  I.  J) +0.  9+DELIN2* 

StlUICEt  I— 1,  J— 1,  K)*ZETA<  I.J>+UICEd,J-l.K>* 

&<-ZETA<  I.  J)+ZETA<  1  +  1.  J)  > -UICEC 1+1.  J-i.  K)#ZETA(  1  +  1.  J> 

&+UICEC  I  — 1.  J,  K)*<  -ZETAC I .  J+l  >+ZETA(  I.  J)  > 

Jc-UICEC  1-1.  J+l.  K)*ZETA(  I.  J+l)  ) 

FX3-0.  9* ( FXETA ( 1 >  +FX ZETA C 1 ) +FXETA ( 2 ) +FXETA < 3 ) +FX  ZETAC  4 ) -FXETA ( 4 ) > 
FX3-FX3+2.  0+THETA 
FXCP*AMASS(I,  J)*THETA* 

&  (UICEC ( I . J)«(UICE( 1  +  1,  J,  1 )-UICE< I  —  1 .  J.  1 >  > 

*<+VICEC(I.  J>*<UICE<I.  J+l.  1  >-UICE(  I.  J-l,  1)))*0.  9*  DEL  IN 
FX3-FX3+FXCP 

FY3-0.  9* <  FYETA ( 1 >  +FYETA <  2  >  +FYZETA <  2 ) +FYZETA ( 3  > -FYETAC3 ) +FYETA <  4 ) > 
FY3*FY3*2.  0*THETA 
FYCP*AMASS(I.  J)*THETA* 

«<  (UICEC  (1<J)*(VICE(I+1>J,  1  )-VICE(  1-1.  J,  1>> 

&+VICEC ( I,  J>*(VICE< I. J+l,  1 >-VICE( I.  J-l.  1 ) ) )*0.  9+DELIN 
FY3-FY3+FYCP 

FL11-0.  3+DRAGAC I,  J)*C0EFI ( I,  J) 

FL1 1“FL1 1*2.  0+THETA 
Fll-CFXMCI.  J)+FX3)*C0EFI < I, J> 

F22*(FYMd,  J)+FY3)*C0EFI < I,  J> 

909  CONTINUE 

FL1 1S=1.  0+FL1 1**2 
FL 1 1 S I  =  1 .  0/FL11S 

UICOR-C  (Fll+FLll#F22)*FLllBI)*UVMd.  J> 

VIC0R=( (F22-FL1 1*F1 1 )*FL1 1SI >*UVM(I.  J) 

OUICE(I,  J,  1 >*0UICE( I,  J,  1 ) +WFA* (UICOR-UICE ( I,  J.  1 >  > 

OVICEd,  J.  1  )*QV1CE(  I.  J,  1  >+WFA*(VICQR-VICE(  I,  J,  1>) 

*  OUICEd,  J.  1  )*OUICE(I,  J,  1  >  +  l.  9*  ( UICOR-UICE  ( I,  J.  1  > ) 

*  OVICEd.  J.  1  > “OVICEd.  J.  1  »  +  l .  9*(VIC0R-VICE<  I.  J.  1>> 

UICEd.  J.  1  >*OUICE(  I.  J.  1> 

VICEd.  J.  1)»0VICE(I.  J,  l> 

*  UICEd.  J.  1 ) “UICEC  I.  J.  1)+1.  9*<UIC0R-UICE(  I.  J.  1>> 

*  VICEd.  J.  1>-VICE(I.  J.  l)  +  l.  3*(  VICOR-VICE<  I,  J,  1>) 

103  CONTINUE 
ICOUNT=ICOUNT+l 
IFdCOUNT.  GT.  BOO)  00  TO  201 
IF < ICOUNT.  LT.  100)  00  TO  202 
WFA*1.  0 

202  CONTINUE 

*  NOW  CHECK  MAX  ERROR 

*  SI 1*0  0 

*  S22*0.  0 

*  FORM  ERROR  MATRIX 

DO  104  J*l. NY 
DO  104  1*1. NX 

UERR < I .  J)*UICE( 1,  J.  1 >-UICE< I,  J, 3) 

VERRd,  J)*VICE(I.  J,  1 ) -VICE (I,  J,  3) 

104  CONTINUE 

*  NOW  FIND  ERROR 

CALL  XMAXM(UERRV) SI ) 

CALL  XMAXMCVERRV.  S2) 

*  81*ABS(UERRV( IU+1 ) ) 

*  82*ABS(VERRV( IV+1 ) ) 

81-DMAXKS1.  S2> 

IF(S1.  LT.  ERROR)  00  TO  200 


25 


11140  00  TO  100 

11150  201  CONTINUE 

11160  PRINT  11 

11170  11  FORMAT <1X.  'NO  CONVERGENCE  AFTER  BOO  ITERATIONS') 

11180  •  NOW  END 
11190  200  CONTINUE 

11200  *  PRINT  1,  I COUNT 
11210  PRINT  12, SI 

11220  PRINT  1,  1 COUNT 

11230  12  FORMATUX,  'MAX  ERROR  AND  U  AND  V  POWER  ',  3012  5) 

11240  1  FORMATUX,  'NO  OF  ITERATIONS  ARE',  1012.  5) 

11250  RETURN 

11260  END 

11270  SUBROUTINE  FELLIP (UICE. VICE, ETA. F, I , J, K ) 

11280  *  SPACER 

11290  IMPLICIT  DOUBLE  PRECISION  (A-H. 0-Z> 

11300  PARAMETER  NX=5, NY-6.  N3-(NX+1 >*(NY+1  ),  NX1-NX+1 

11310  it,  NYl-NY+l.NXMl-NX-l.  NYM1-NY-1 

11320  DIMENSION  UICE (NX.  NY,  3) ,  VICE (NX,  NY,  3),  ETA ( NX l .  NY1 ), F(4) 

1 1330  COMMON/STEP /OTAT,  OTAX,  OTAY 

11340  DELTAX-OTAX 

11350  S1-.  5/(DELTAX**2> 

11360  F ( 1 > -SI* ( UICE ( I+l. J,K)*<ETA<I+1. J+l >+ETA( I  +  l, J) ) 

11370  &-UICE( I,  J.  K)*(ETA( I  + 1 »  J+ 1  > +ET A ( I , J >  +ET A ( I  +  l. J)+ETA( I.  J+t ) > 

11380  &+UICEI  1-1, J. K)*(ETA( I.  J+l  >+ETAU.  J)  )  ) 

11390  F(2>=S1*(UICE(I,  J+l ,  K ) * < ETA( I  +  l.  J+1>+ETA(I,  J+l>  > 

11400  &-U1CEU.  J,  K)*(ETA(1  +  1,  J+l >+ETA(I,  J)+ETA(  I  +  l . J>+ETA< 1 .  J+l  >  ) 

11410  S<+UICE( 1,  J-l.  K)*(ETA( I,  J)+ETA( I  +  l, J) ) ) 

11420  F(3)-S1*(VICE( I  — 1 , J-l,  K)#ETA( I,  J)+VICE(I, J-l. K>*(-ETA< I. J) 

11430  &+ETA(  I  +  l,  J)  ) -VICE ( I+l,  J-l,  K ) *ETA(  I+l.  J)+VICE<  1-1.  J,  K)*(-ETA(  I,  J+l ) 

11440  &+ETAU,  J>  )  +VICE(I.  J.  K)*(-ETA(  I.  Jl-ETAU  +  l.  J+l >+ETA( 1+1.  J) 

11450  (k+ETA(  I. J+l ) ) ) 

11460  F(3)=F(3)+S1*(VICE(  I  +  l.  J,  K)+(-ETA(  I  +  l ,  J>+ETA(  I  +  l ,  J+l  >  > 

11470  Sc-VICE(  1-1.  J+l.  K)*ETA<  I,  J+l  > 

11480  &+VICE( I.  J+l.  K)#(-ETA( I  +  l. J+1)+ETA(I, J+l >  ) 

11490  *.+VICE(  I  +  l.  J+l.  K)*ETA(  I  +  l,  J+l )  ) 

11500  F ( 4 )— SI* (VICE ( I— 1,  J-l,  K)*ETA( I.  J)+VICE( I.  J-l.  K)*(-ETA( I  +  l.  J> 

11510  *c+ETA( I. J) )-VICE( I  +  l.  J-l.  K>*ETA( I+l. J>+VICE(I-1. J. K) *(ETA( I , J+l ) 

11520  6-ETA ( I . J) )+VICE( I.  J,  K)*(ETA( I.  J+l )+ETA( I  +  l. J)-ETA( I.  J)-ETA( I  +  l, 

11530  «<J+1>>) 

11540  F(4)»F(4)+S1*(VICE(I  +  1.  J.  K ) *<ETA< I  +  l , J>-ETA( I+l .  J+l >  > 

11550  Si-VICE (  1-1,  J+l. K)*ETA(I. J+l) 

11560  «<+VICE(  I.  J+l. K)*(ETA<  I  +  l,  J+l  >-ETA(  I,  J+l  )  ) 

11570  6+VICEU  +  l.  J+l.  K)*ETA(  I  +  l,  J+l ) ) 

11580  F(3)-F(3)  *.  5 

11590  F(4)-F(4)*.  5 

11600  RETURN 

11610  END 

11620  SUBROUTINE  FELLD1 (UICE, VICE. ETA,  F,  K ) 

11630  •  SPACER 

11640  IMPLICIT  DOUBLE  PRECISION  (A-H, O-Z) 

11650  PARAMETER  NX-5. NY-6.  N3-(NX  +  1 )*<NY+1 >.  NX 1-NX+l 

11660  6.  NY1-NY+1,  NXM1-NX-1.  NYM1-NY-1 

11670  DIMENSION  UICE(NX,  NY.  3),  VICEINX, NY, 3) , ETA(NX1 , NY1 > .  F ( 4.  NX. NY > 

11680  *  COMMON/STEP/DEL TAT,  DEL TAX. DELTAY 

11690  COMMON/DINV/Sl 

11700  DO  101  J-2.  NYM1 

11710  DO  101  1-2, NXM1 

11720  F( 1,  I.  J>-S1*(UICE( I  +  l.  J, K)*(ETA< I  +  l,  J+l )+ETA( I  +  l, J> ) ) 

11730  F(2. I. J>-S1*(UICE( I.  J+l.  K)*(ETA< I  +  l,  J+l J+ETAU.  J+t> ) ) 

11740  F<3. I. J)-S1*(VICE( I,  J,  K)*(-ETA(I.  J)-ETA( I  +  i.  J+t )+ETA< I  +  l,  J> 

11750  li+ETAU.  J+l)) 

11760  fc+VICEI I  +  l.  J. K)*( -ETA( I  +  l ,  J)+ETA( I  +  l,  J+l ) ) 

11770  lc+VICE(  I,  J+l,  K)*(-ETA(  I  +  l ,  J+ 1 ) +ETA<  I ,  J+l  >  > 

11780  l«+VICE<  I  +  l.  J+l. K)  *ETA(  I  +  l.  J+l) ) 

11790  F<4.  I.  J)-Sl*<VICE( I. J, K)*(ETA( 1 , J+l ) +ETA( I  +  l. J)-ETA( I, J)-ETA( I  +  l . 


11800  1<J+1 ) >+VICE< I  +  l.  J.  K)#<ETA<  I  +  l,  J»-ETA< I  +  l.  J+l ) ) 

11810  l«+VICE<  I<  J+l  -  K)*<ETA<  X  +  l-  J+l >-ETA<I.  J+l  )  ) 

11820  li+VICE ( 1  +  1 > J+l.  K)*ETA< I  +  l.  J+l >  > 

11830  F<3.  I.  J>-0  9*F<3.  I.  J) 

11840  F<4.  I. J)-0  9*F<4.  I.  J> 

11890  101  CONTINUE 

11860  RETURN 

11870  END 

11880  SUBROUTINE  F£LLD<UICE,  VICE.  ETA.  F.  I.  J.  K. FETA) 

11890  *  SPACER 

11900  IMPLICIT  DOUBLE  PRECISION  <A~H. 0-Z> 

11910  PARAMETER  NX-9.  NY-6.  N3-<NX+t )* <NY+1 ), NX1-NX+1 

11920  t<.  NV1-NY+1 ,  NXM1-NX-1,  NVM1-NY-1 

11930  DIMENSION  UICE< NX,  NY.  3 > .  VICE<NX.  NY,  3 ) .  ETA<NXl,  NY1  > ,  F < 4> 

11940  fc,  FETA <4, NX, NY) 

11990  COMMON/ D I NV/S1 

11960  F< 1 >-FETA< 1. I, J)+S1«< 

11970  fc+UICE< I - 1 .  J.  K)«<ETA< I.  J+l )+ETA( I, J) ) ) 

11980  F<2)«FETA<2, I, J)+S1#< 

11990  &+UlCE< I,  J— 1.  K)*<ETA< I.  J)+ETA< I  +  l.  J) >  > 

12000  F<3)-FETA<3, I. J)+0  9*S1*<VICE< I - 1 .  J-l. K)*ETA< I. J)+VICE< I. J-l.K)* 

12010  «<(-ETA(I,J> 

12020  *<+£TA<  I  +  l.  J) >-VICE< I+l.  J-l,  K)*ETA(  I  +  l,  J)+VICE<  1-1.  J.  K)*<-ETA<  I.  J+l  > 

12030  1<+ETA<  I,  J) ) 

12040  &-VICE< 1-1.  J+l.  K)*ETA<  I.  J+t )  > 

12050  F<4)*FETA<4.  I. J)+S1*0  5* 

12060  *<< VICE <1-1.  J-l.  K)*ETA<  I.  J)+VICE<  I,  J-l ,  K >#< -ETA<  I  +  l.  J) 

12070  &+ETA< I,  J) > -VICE < 1  +  1,  J-l.  K)*ETA< 1*1 , J)+VICE< 1-1. J, K)*CETA( I, J+l > 

12080  *c-£TA(I.J>) 

12090  gr-VICE< 1-1 . J+l , K  > *ETA< 1 ,  J+l>> 

12100  RETURN 

12110  END 

12120  SUBROUTINE  XMAXM(U£RRV.  X) 

12130  *  SPACER 

12140  IMPLICIT  DOUBLE  PRECISION  <A-H, 0~Z> 

12150  PARAMETER  NX-5.  NY-6.  N3»<NX+1 )*<NY+1 ),  NX1-NX+1 

12160  &, NY1-NY+1 , NXM1-NX-1.  NYM1— NY-1,  N4*NX*NY 

12170  DIMENSION  UERRV<N4) 

12180  X-0  0 

12190  DO  100  I-1.N4 

12200  X-DMAX 1<DABS<UERRV<  I >).  X) 

12210  100  CONTINUE  *  ' 

12220  RETURN 

12230  END 

12240  SUBROUTINE  XSUM < HEFF 1 . S 1 ) 

12250  *  PROGRAM  SUMS  UP  VECTOR 

12260  IMPLICIT  DOUBLE  PRECISION  <A-H, 0-Z) 

12270  PARAMETER  NX =5,  NY-6. N3- (NX+l >  *<NY+1 > , NX 1-NX  +  l 

12280  NYl-NY+l, NXM1-NX-1,  NYM1-NY-1 

12290  DIMENSION  HEFF1 <N3> 

12300  Sl-O.  0 

12310  DO  100  1*1. N3 

12320  S1-S1+HEFF1 < I ) 

12330  100  CONTINUE 

12340  RETURN 

12350  END 


APPENDIX  B:  TWENTY-ONE-DAY  TEST  RUN  OF  MODEL 

This  21 -day  test  run  is  obtained  by  compiling  and  executing  the  code  as  listed  in  Appendix  A.  The 
only  input  data  needed  are  contained  in  a  file  BNDATA  which  is  read  by  subroutine  BNDRY  (lines 
1020-1050).  This  BNDATA  file  consists  of  the  following  formatted  data  which  define  the  boundary 
conditions  shown  in  Figure  2: 

000000 11100111 00 11100111 000000 
0000000 1111001111001111 00 1111001111 0000000 
0000000 1 1 1 000 1 1 11001 1 1 10011 11001111 0000000 


INTEGER  VARIABLE  (K2)  SET,  NEVER  USED,  LINES  250  AND  300 


DOUBLE  PRECISION  VARIABLE 
DOUBLE  PRECISION  VARIABLE 
DOUBLE  PRECISION  VARIABLE 
DOUBLE  PRECISION  VARIABLE 
DOUBLE  PRECISION  VARIABLE 


(DWAT)  SET,  NEVER  USED.  LINE  2660 
(DAIR)  SET,  NEVER  USED.  LINE  2670 
(CRAV>  SET.  NEVER  USED,  LINE  2750 
(EPSIL)  SET,  NEVER  USED,  LINE  3210. 
<ETAC )  SET,  NEVER  USED,  LINE  3240 


DOUBLE  PRECISION  VARIABLE  (ONE)  SET.  NEVER  USED,  LINE  3890 
DOUBLE  PRECISION  VARIABLE  (ZERO)  SET,  NEVER  USED,  LINE  3900 


INTEGER  VARIABLE  <KGR0>  SET,  NEVER  USED,  LINE  6200 
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1  OOOOOOOOOOOOO 
I  OOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 
1  ooooooooooooo 
1  ooooooooooooo 


OOOOOOOOOOOOOO 

1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooouooooo 
i  ooooooooooooo 

OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
1 .  OOOOOOOOOOOOO 

1  ooooooooooooo 

1  OOOOOOOUOOOOO 
l  UUUUUUOOOOOOO 

1  OOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 

1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 


OOOOOOOOOOOOOO 

1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 
I  ooooooooooooo 

OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 
1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 


OOOOOOOOOOOOOO 
1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 
i  ooooooooooooo 


OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
' OOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 
1  ooooooooooooo 


WESUM  IS  19  0000 

HESUM  is  20  OOOO 

3  2000000000000 
3  3000000000000 
3  4000000000000 
3  5000000000000 
3  6000000000000 
3  7000000000000 


UVSUM  IS  12 

UVSUtt  IS  12 

3  3000000000000 
3  4000000000000 
3  5000000000000 
3  6000000000000 
3  7000000000000 
3  8000000000000 


00000 

00000 

3. 4000000000000 
3  5000000000000 
3  6000000000000 
3  7000000000000 
3  8000000000000 
3  9000000000000 


3  5000000000000 
3  6000000000000 
3  7000000000000 
3  8000000000000 

3  9000000000000 

4  ooooooooooooo 


3  6000000000000 
3  7000000000000 
3  8000000000000 

3  9000000000000 

4  ooooooooooooo 
A  1 oooooooooooo 


3  6000000000000 
3  4000000000000 

j  2000000000000 

3  ooooooooooooo 
-2  8000000000000 
■2  6000000000000 


-3  4000000000000 
-3  2000000000000 

-3  ooooooooooooo 
-2  8000000000000 
-2  6000000000000 
-2  4000000000000 


-3  2000000000000 
-3  ooooooooooooo 
-2  8000000000000 
-2  6000000000000 
-2  4000000000000 
-2  2000000000000 


-3  ooooooooooooo 
-2  8000000000000 
-2  6000000000000 
-2  4000000000000 
-2  2000000000000 
-2  ooooooooooooo 


-2  8000000000000 
-S  6000000000000 
-2  4000000000000 
-2  2000000000000 
-2  ooooooooooooo 
-1  8000000000000 


I OOOGOOOOOOOOOE-O 1 

i looooooooooooe-oi 

I 2000000000000E-01 
1 3000000000000E -O I 
I 4000000000000E-01 
I 5000000000000E-01 


90000000000000E-02 
l OOOOOOOOOOOOOE-O 1 
l 1 OOOOOOOOOOOOE -0 1 
1 2000000000000E -0 1 
1 3000000000000E-0 l 
1 4000000000000E-0 1 


80000000000000E-02 
900000U0000000E-02 
1 OOOOOOOOOOOOOE-O 1 
I  1 OQOOOOQOOOOOE-O I 
I2000000000000E-01 
1 30000 000 000 00E-01 


7 OOOOOOOOOOOOOE -02 
80000000000000E-02 
90000000000000E -02 
t  OOOOOOOOOOOOOE - 0 1 
1 1 OOOOOOOOOOOOE -01 
1 2000000000000E - 0 1 


60000000000000E  -  <j2 

70000000000000*  OJ* 
80000000000000*  -02 
90000000000000F -02 
1 OOOOOOOOOOOOOE -0 1 
1 1 OOOOOOOOOOOOE - t) 1 
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1 


1  1399305555S6E-05 

2  225694444444E-06 

3  532407407407E-07 

4  42824074074  IE-07 

5  35879629 6296E-07 

6  3 1 2500000000E-07 

7  243055555556E-07 

8  1 62037037037E  -07 

9  104166666667E-07 

10  347222222222E-0B 

1 1  000000000000 

12  -  231 48 1481 48 IE— 08 

13  -  347222222222E-08 

14  -  462962962963E-08 

15  -  578703703704E-08 

16  -  578703703704E-08 

17  -  694444444444E-08 

18  -  694444444444E-08 

19  -  694444444444E-08 

20  -  694444444444E-08 

21  -  694444 444444E -00 


9 1 984855243 125E -02 
4851 4339437 489E -0 1 
491001 200387 5 1 E-0 1 
496 135311 69529E-0 1 
5005 1 5 1 353279BE-0 1 
1 1 202822487627E-0 1 


16909 1 33238402E-0 1 
44866731 19B317E-02 
374323493 1 08 1 8E -02 
307489934697 16E-02 
24843979837337E-02 
.  1 467603947 1344E-01 


1 5793329942552E-01 
.  32 1 37460399924E-02 
.  253490056S8436E-02 
19335467202533E-02 
141201 33S32S38E-02 
1 3893 1 39793508E-0 1 


.  537 1 9427407723E-01 
4 1 095469 1 4 1 595E-0 1 
.  1 3850900 150383E-02 
.  8532064 5905588E-03 
.  40264520065 144E-03 
131 7472347 1 373E-0 l 


2 598535324430 2£ -01 
.  52528457534 129E -01 
7893 1 630880927E -0 1 
7869022 1 589046E-0 1 
78531 100243467E-01 
.  52167493728903E-01 


-  312701 4647 5406E-0 1 - 
5624934887490 7E-02 
84022932 1 277.0 1 E-02 
1 1 1 33006 1 49244E-0 1 
1 38267666 1 20 1 9E-0 1 
4781 B393763294E-0 1 
MAX  ERROR  AND  U  AND  V 
NO  OF  ITERATIONS  ARE 


61104461 2497 1 5E-0 1 - 
97804294854293E-02 
1 27 1 352533670 1 E-0 1 
1 5609S95065236E-0 1 
1 S478933637849E— 01 
842004 1 5473009E— 0 1 
POWER  31432E— 05 

11  000 


59131 530 18387 1E-01- 
1 147 1 2363 l 5805E-01 - 
.  1435081 5038782E-01 
1720357771 1480E-01 
.  200403424 50 109E-01 
85959404 535894E— 01 


250 1 39645S267BE-0 1 
20727131 162579E-01- 
1595231 1038522E -01 
1 B77262592 1 528E-0 1 
21 S88086253953E-01 
8771 62B7349946E-0 l 


927624277538 1 4E-02 
2054196541 1089E-01 
1 5776677375776E-0 1 
.  1 8356824 5284 1 OE-0 1 
.  20943549621 534E-01 
5531 1079B04659E— 01 


00000000000000 

00000000000000 

00000000000000 

00000000000000 

oooooooooooooo 

00000000000000 


oooooooooooooo 
55 1 23024654 1 73E-03 
3939339 1 967 1 57E-03 
1 9239739927943E-03 
-.  14394697336957E-05 
OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
.  B2694462070983E-03 
.  499955 1945872 1 E— 03 
2758 1 6 1 096O687E-03 
.  1261 1643083897E-03 
.  OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
1 2280766662 1 05E-02 
296655774727B5E— 03 
.  191 9484 18S3472E-03 
.  19914970340528E-03 
OOOOOOOOOOOOOO 


.  OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 

oooooooooooooo 
.  oooooooooooooo 
.  oooooooooooooo 
oooooooooooooo 


oooooooooooooo 
oooooooooooooo 
oooooooooooooo 
oooooooooooooo 
oooooooooooooo 
.  oooooooooooooo 


.  oooooooooooooo 

.  6446B343209B46E-03 
954 1 75029899 1 4E-03 
.  1052821 1003075E-02 
.  790484752521 98E -03 
.  OOOOOOOOOOOOOO 


.  OOOOOOOOOOOOOO 
.  70332404696872E-03-. 
.  1 18751 37686298E-02 
.  1 4008863763407E-02  . 
.  10379 124877481 E-02 
.  OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
93075828B95523E— 04 
842443984 1 7 1 97E-03 
1 17769 13525051 E-02 
93501 620599902E-03 
OOOOOOOOOOOOOO 


.  OOOOOOOOOOOOOO 
.  OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 

.  oooooooooooooo 
oooooooooooooo 
oooooooooooooo 


MAX  ERROR  AND  U  AND  V  POWER  .  23423E-05 
NO  OF  ITERATIONS  ARE  11.000 
MAX  ERROR  AND  U  AND  V  POWER  40834E-05 
NO  OF  ITERATIONS  ARE  8.0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  1.  OOOOOOOOOOO  3.  29669230773 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  1.00000000000  62.6692788514 
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OUTFLOW  AND  NET  OUTFLOW  ARE  .  9147826E-04  9147826E-04 


OOOOOOOOOOOOOO 
.  OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
.  OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
41057247330924E-04 
28 1 389292004 7BE-04 
1 243 1 366 1 069 1 0E-04 
-  17442721 54909 IE-05 
.  OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
36846033029838E-04 
33823664003268E-04 
1 864687834 1 670E-04 
.  93033439 1676 7 IE-03 
.  OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
8 1 4680 38298827E-04 
1848241 101 1226E-04 
.  1308343 173 343BE-04 
.  145331 202 1 3 194E-04 
OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 


,  OOOOOOOOOOOOOO  OOOOOOOOOOOOOO  OOOOOOOOOOOOOO 

.  34421 302682481E-04  .  61 76622343294 IE-04  .  1 3906653062933E-03 
82072462036398E-04  .  103024036 17843E-03  .  74681 43274743BE-04 
.  90320023744 148E-04  .  120 790278723 19E-03  .  10016494703124E-03 
.  66972 1701 70303E -04  .  B9820544754459E-04  797013981 181 3PE-04 

.  OOOOOOOOOOOOOO  OOOOOOOOOOOOOO  .  OOOOOOOOOOOOOO 


.  OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
.  OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 


.  OOOOOOOOOOOOOO 
3.  2982791312133 
3.  2982775716448 
3.  2983322686140 
3.  2984023133110 
3.  2984662040978 


OOOOOOOOOOOOOO 
.  99996700313530 
.  99996633006558 
99998312147937 
1 .  0000000000000 
1 .  OOOOOOOOOOOOO 


.  OOOOOOOOOOOOOO 
3.  2982373791387 
3.  2982826576149 
3.  2983470141 199 
3.  29843014931 17 
3.  2985539676780 


OOOOOOOOOOOOOO 
.  99718959907436 
.  99996807280986 
99998739428303 
1 .  OOOOOOOOOOOOO 
1 .  OOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
3  2982877200233 
3  2982446106205 
3  2983647286222 
3  2984468663240 
3  2983730780206 


OOOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 
99442693180934 
99999296768207 
1  OOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
.  OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
. OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
3  2984181030719 
3  2983948719993 
3  2984427139725 
3.  2984933018081 


OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 
99723731117466 
1  OOOOOOOOOOOOO 
1  OOOOOOOOOOOOO 


OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
4120879120879.  1 
41 20879 120B79  1 
4120879120879  1 
4120879120879.  1 
4120879120879  1 


OOOOOOOOOOOOOO 
4120879120879  1 
4120079120879  1 
4120879120879  1 
4120879120879  1 
4120879120879.  1 


OOOOOOOOOOOOOO 
4120879120879  1 
3910418381531.  3 
4120879120879.  1 
4120879120879.  1 
3630562648686.  7 


OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
4120879120879.  1 
4120879120879  1 
4120879120879  1 
4120879120879.  1 


OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 


OOOOOOOOOOOOOO 
1 6846 1 33846 1 54E -02 
1 6846 1 33846 1  54E -02 
.  16846 1 33846 1 34E-02 
1 6846 1 33846 1 34E-02 
1 6846 1 53846 1 54E-02 


.  OOOOOOOOOOOOOO 
.  168461  538461  34E-02 
1 68461 53846 1 54E-02 
1 6846 1 33846 1 34E-02 
.  168461 33846 134E-02 
.  168461 33846 154E-02 


.  OOOOOOOOOOOOOO 
1 6846 1 33846 1 34E -02 
1 6846 1 33846 1 34E -02 
1 6846 1 33846 1 34E-02 
1 6846 1 53846 1 34E-02 
1 6846 1 33846 1 34E-02 


.  OOOOOOOOOOOOOO 
OOOOOOOOOOOOOO 
1 6846 1 53846 1 34E-02 
.  168461 33846 134E-02 
.  168461 33846 134E-02 
.  1 68461 33846 134E -02 


OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 

OOOOOOOOOOOOOO 


.  OOOOOOOOOOOOOO 
-3941.  0269320377 
-3338  6697984968 
-3743.  3721734968 
-7943.  2686226044 
-10397  702303097 


OOOOOOOOOOOOOO 
-2998  2341796073 
-4740.  3630366613 
-7692  1710279809 
-8339  3741244123 
-13449  670173992 


.  OOOOOOOOOOOOOO 
-4613  0152603708 
-2960  4273986860 
-8373  7826942921 
-9082.  1036741092 
-13633  964319437 


.  . OOOOOOOOOOOOOO 
-6808461394898  3 
-11643  622234333 
-9389.  7083890993 
-10038.  926290415 
-11813.  303386032 


1 3993033333336E-03  .  22369444444444E-06 
3 1 230000000000E-07  24303333333336E-07 


.  3324074074074 IE-07 
.  1 6203703703704E-07 


42824074074074E— 07  .  33B79629629630E-07 
104 1 6666666667E-07  .  34722222222222E-08 


W1K  UWWWVWWVt  V  r  ,  wwr  —  '  •  -  —  ’  -  —  —  — — —  —  -  —  - 

OOOOOOOOOOOOOO  -.  23 1481 48 1 48 1 4BE-08- .  34722222222222E-08-.  46296296296296E-08-.  37870370370370E-08 
37870370370370E-08-.  69444444444444E-08-.  69444444444444E-08-.  69444444444444E-08-.  69444444444444E-08 
69444 444444444E-08 


SQUARE  VELOCITY,  80  VELOCITY  DIFFERENCE, MAX  CHANOE 

9661328361 19E-07  1 1943300 1931 E-04  .  128009609762E-02 
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TIME  STEP  AND  TOTAL  THICKNESS  ARE  1  OOOOOOOOOOO  62  6692788314 

OUTFLOW  AND  NET  OUTFLOW  ARE  9147826E-04  9147826E-04 

GROWTH  AND  NET  GROWTH  ARE  3200769E-01  .  3200769E-01 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  0000000  0000000  18  98872 

MAX  ERROR  AND  U  AND  V  POWER  41442E-03 

NO  OF  ITERATIONS  ARE  2  0000 

MAX  ERROR  AND  U  AND  V  POWER  24269E-05 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  2.  OOOOOOOOOOO  62.  7024471339 

OUTFLOW  AND  NET  OUTFLOW  ARE  8751464E-04  .  1789929E-03 

GROWTH  AND  NET  GROWTH  ARE  .  33255B0E-01  .  6326349E-01 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  .  1363B57E-02  1363857E-02  18  99109 

MAX  ERROR  AND  U  AND  V  POWER  2B397E-05 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  13899E-05 

NO  OF  ITERATIONS  ARE  1.0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  3.  OOOOOOOOOOO  62.  7353028718 

OUTFLOW  AND  NET  OUTFLOW  ARE  . 8244857E-04  . 2614415E-03 

GROWTH  AND  NET  GROWTH  ARE  .  32938 19E-01  .  9820168E-01 

OPEN  WATER  CROWTH  AND  NET  GROWTH  ETC  .  1077774E-02  2441631E-02  18  99294 

MAX  ERROR  AND  U  AND  V  POWER  .  96390E-06 

NO  OF  ITERATIONS  ARE  1.0000 

MAX  ERROR  AND  U  AND  V  POWER  88624E-06 

NO  OF  ITERATIONS  ARE  1.0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  4.  OOOOOOOOOOO  62.  7679033279 

OUTFLOW  AND  NET  OUTFLOW  ARE  .  7928826E-04  .  3407297E-03 

GROWTH  AND  NET  GROWTH  ARE  3267974E-01  .  1308814 

OPEN  WATER  GROWTH  AND  NET  CROWTH  ETC  S336313E-03  3293283E-02  18.99438 

MAX  ERROR  AND  U  AND  V  POWER  . 66133E-06 

NO  OF  ITERATIONS  ARE  1.0000 

MAX  ERROR  AND  U  AND  V  POWER  67121E-06 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  5.  OOOOOOOOOOO  62  8002938478 

OUTFLOW  AND  NET  OUTFLOW  ARE  7855221E-04  4192819E--03 

GROWTH  AND  NET  CROWTH  ARE  3246907E-01  .  1633505 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  .  6795474E-03  3974830E-02  1 8  99349 

MAX  ERROR  AND  U  AND  V  POWER  46274E-06 

NO  OF  ITERATIONS  ARE  1.0000 

MAX  ERROR  AND  U  AND  V  POWER  45260E-06 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  6  OOOOOOOOOOO  62  0325112472 

OUTFLOW  AND  NET  OUTFLOW  ARE  792U90E-04  49B4938E-03 

GROWTH  AND  NET  GROWTH  ARE  .  3229661E-01  .  1956471 

OPEN  WATER  GROWTH  AND  NET  CROWTH  ETC  3434454E-03  .  4520276E-02  18  99634 

MAX  ERROR  AND  U  AND  U  POWER  36399E-06 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  .  31219E-06 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  7  OOOOOOOOOOO  62  B645B45271 

OUTFLOW  AND  NET  OUTFLOW  ARE  8060446E-04  .  S790983E-03 

GROWTH  AND  NET  GROWTH  ARE  321S3B8E-01  .2278010 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  4424361E-03  4962732E-02  18  99699 

MAX  ERROR  AND  U  AND  V  POWER  28442E-06 

NO  OF  ITERATIONS  ARE  1.0000 

MAX  ERROR  AND  U  AND  V  POWER  21306E-06 

NO  OF  ITERATIONS  ARE  1.  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  8.  OOOOOOOOOOO  62  8963363328 

OUTFLOW  AND  NET  OUTFLOW  ARE  .  8238654E-04  .  6614848E-03 

OROWTH  AND  NET  GROWTH  ARE  .  32034 19E-01  2598332 

OPEN  WATER  OROWTH  AND  NET  OROWTH  ETC  .  363SS32E-03  3326287E-02  18  99748 

MAX  ERROR  AND  U  AND  V  POWER  22296E-06 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  .  14870E-06 

NO  OF  ITERATIONS  ARE  1.0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  9.  OOOOOOOOOOO  62.  9283834123 

OUTFLOW  AND  NET  OUTFLOW  ARE  8433947E-04  .  743B243E-03 

GROWTH  AND  NET  OROWTH  ARE  3193342E-01  .2917686 

OPEN  WATER  OROWTH  AND  NET  OROWTH  ETC  .  3044336E-03  .  3630721E-02  18  99783 

MAX  ERROR  AND  U  AND  V  POWER  .  1761BE-06 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  .  10368E-06 

NO  OF  ITERATIONS  ARE  1.0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  10  OOOOOOOOOO  62.  9601434288 

OUTFLOW  AND  NET  OUTFLOW  ARE  8631673E-04  .  B321410E-03 

OROWTH  AND  NET  OROWTH  ARE  3184633E-01  .  3236149 

OPEN  WATER  OROWTH  AND  NET  OROWTH  ETC  .  2396003E-03  3890321E-02  18  99813 

MAX  ERROR  AND  U  AND  V  POWER  .  14327E-06 


18  98872 


37  97980  | 

I 

i 

j 

56.  97274  j 


73  96712 


94  96261 


113  9590 


132  9339 


131  9334 


170  9313 


189  9494 
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NO  OF  ITERATIONS  ARE  1  OOOO 

MAX  ERROR  AND  U  AND  V  POWER  74802E-07 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  11  0000000000  62  9918267462 

OUTFLOW  AND  NET  OUTFLOW  ARE  8824076E-04  9203818E-03 

GROWTH  AND  NET  GROWTH  ARE  3176956E-01  3553845 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  2255131E-03  6U5834E-02  18  99835 

MAX  ERROR  AND  U  AND  V  POWER  11812E-06 
NO  OF  ITERATIONS  ARE  1  0000 
MAX  ERROR  AND  U  AND  V  POWER  60845E-07 
NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  12  0000000000  63  0234373316 

OUTFLOW  AND  NET  OUTFLOW  ARE  9006798E-04  1010450E-02 

GROWTH  AND  NET  GROWTH  ARE  3170065E-01  3870851 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  1996125E-03  6315447E-02  18  99851 

MAX  ERROR  AND  U  AND  V  POWER  99636E-07 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  53976E-07 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  13  0000000000  63.0549833090 

OUTFLOW  AND  NET  OUTFLOW  ARE  9177519E-04  1102225E-02 

GROWTH  AND  NET  GROWTH  ARE  3163775E-01  4187229 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  .  1799415E-03  6495388E-02  18  99864 

MAX  ERROR  AND  U  AND  V  POWER  85949E-07 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  49112E-07 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  14  0000000000  63  0064693874 

OUTFLOW  AND  NET  OUTFLOW  ARE  933532BE-04  1195578E-02 

GROWTH  AND  NET  GROWTH  ARE  3157943E-01  4503023 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  1650067E-03  6660395E-02  18  99873 

MAX  ERROR  AND  U  AND  V  POWER  80681E-07 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  45665E-07 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  1 5  0000000000  63  1178991909 

OUTFLOW  AND  NET  OUTFLOW  ARE  9480144E-04  1290380E-02 

GROWTH  AND  NET  GROWTH  ARE  .  3152460E-01  4818269 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  1 536697E-03  6814064E-02  18  99880 

MAX  ERROR  AND  U  AND  V  POWER  .  76632E-07 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  .  43203E-07 

NO  OF  ITERATIONS  ARE  1.0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  16.  0000000000  63  1492755131 

OUTFLOW  AND  NET  OUTFLOW  ARE  .  9612368E-04  1386503E-02 

GROWTH  AND  NET  GROWTH  ARE  3147245E-01  .  5132994 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  .  1450639E-03  .  695912BE-02  18  99885 

MAX  ERROR  AND  U  AND  V  POWER  .  73288E-07 

NO  OF  ITERATIONS  ARE  1.0000 

MAX  ERROR  AND  U  AND  V  POWER  41433E-07 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  17.  0000000000  63.  1806005122 

OUTFLOW  AND  NET  OUTFLOW  ARE  .  9732675E-04  1483830E-02 

GROWTH  AND  NET  GROWTH  ARE  .  3142233E-01  .  5447217 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  .  1385303E-03  .  709765BE-02  18  99889 

MAX  ERROR  AND  U  AND  V  POWER  .  69765E-07 

NO  OF  ITERATIONS  ARE  1.  0000 

MAX  ERROR  AND  U  AND  V  POWER  .  40019E-07 

NO  OF  ITERATIONS  ARE  1.0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  18.0000000000  63.2118759303 

OUTFLOW  AND  NET  OUTFLOW  ARE  98420 14E-04  .  1 582250E-02 

GROWTH  AND  NET  GROWTH  ARE  .  31373B4E-01  .  5760955 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  .  1336429E-03  .  7231301E-02  18  99892 

MAX  ERROR  AND  U  AND  V  POWER  67617E-07 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  3B332E-07 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  19.  0000000000  63  2431035050 

OUTFLOW  AND  NET  OUTFLOW  ARE  9942083E-04  16B1671E-02 

GROWTH  AND  NET  GROWTH  ARE  3132700E-01  6074225 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  .  1304179E-03  7361719E-02  18.99894 

MAX  ERROR  AND  U  AND  V  POWER  67242E-07 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  30222E-O7 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  20.0000000000  63  2742841446 

OUTFLOW  AND  NET  OUTFLOW  ARE  1003302E-03  1762001E-02 
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265  9436 


284  9424 
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GROWTH  AND  NET  GROWTH  ARE  3128097E-01  6387033 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  1279838E-03  7489703E-02  10  99896  379  9370 

MAX  ERROR  AND  U  AND  V  POWER  66773E-07 

NO  OF  ITERATIONS  ARE  1  0000 

MAX  ERROR  AND  U  AND  V  POWER  36318E-07 

NO  OF  ITERATIONS  ARE  1  0000 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  21  0000000000  3  32938600933 

TIME  STEP  AND  TOTAL  THICKNESS  ARE  21  0000000000  63  3054187372 


OUTFLOW  AND  NET  OUTFLOW  ARE  311530E-03  18B3154E-02 
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ftf  rtf  rtf 


13993039399336E-09  22969444444444E-06  93240740740741E-07  42824074074074E-07  39879629629630E-07 
3 1 290000000000E-07  24309999999996E-07  .  16203703703704E-07  10416666666667E-07  34722222222222F-08 
.  00000000000000  -  23148148148148E-08-  34722222222222E-08-  46296296296296E-08-  97870370370370E-08 

-  37870370370370E-08-  69444444444444E-08-.  69444444444444E-08-  69444444444444E-0B-  69444444444444E-08 

-  6944444 44444 44E-08 


SQUARE  VELOCITY, SQ.  VELOCITY  DIFFERENCE,  MAX  CHANCE 

B93987030642E-0?  .  698437122999E-13  847259876693E-07 


TIME  STEP  AND  TOTAL  THICKNESS  ARE  21  0000000000  63  3034187372 

OUTFLOW  AND  NET  OUTFLOW  ARE  10U330E-03  18831  94E-02 

GROWTH  AND  NET  OROWTH  ARE  .  3123973E-01  .  6699393 

OPEN  WATER  GROWTH  AND  NET  GROWTH  ETC  1263262E-03  7616029E-02  18  99897  398  9360 

14.  277  soc  901  I/O 
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